sn-ENS neuron integration
advanced:
clustering modification
cell composition
all markers
check signature score of some genesets
source("/Shared_win/projects/RNA_normal/analysis.10x.r")load integrated obj
GEX.seur <- readRDS("./sn10x_WYS.sct_anno.rds")
GEX.seur## An object of class Seurat
## 47356 features across 19649 samples within 3 assays
## Active assay: integrated (2397 features, 2397 variable features)
## 2 other assays present: RNA, SCT
## 3 dimensional reductions calculated: pca, tsne, umap#scales::show_col(ggsci::pal_igv("default")(49))color.cnt <- scales::hue_pal()(4)[c(2,1,4,3)]
color.cnt## [1] "#7CAE00" "#F8766D" "#C77CFF" "#00BFC4"head(GEX.seur@meta.data)## orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCACAAGACGAC-1_1 Nb5d.PBS_INF 3257 1801 0.36843721 0.3991403
## AAACCCAGTGGGCTCT-1_1 Nb5d.PBS_INF 1511 997 0.66181337 0.4632694
## AAACCCAGTTTGTTCT-1_1 Nb5d.PBS_INF 2855 1577 0.98073555 0.3152364
## AAACCCATCCTAGCCT-1_1 Nb5d.PBS_INF 2433 1451 0.08220304 0.3699137
## AAACCCATCGAAACAA-1_1 Nb5d.PBS_INF 3129 1656 0.12783637 0.4474273
## AAACCCATCGGTCAGC-1_1 Nb5d.PBS_INF 2201 1294 0.22716947 0.2271695
## S.Score G2M.Score Phase cnt rep newAnno
## AAACCCACAAGACGAC-1_1 0.011590405 -0.0004169865 S Nb5d.INF rep4 EMN3
## AAACCCAGTGGGCTCT-1_1 -0.024203070 0.0012414826 G2M Nb5d.PBS rep4 IPAN1
## AAACCCAGTTTGTTCT-1_1 -0.013980476 0.0039329410 G2M Nb5d.INF rep1 EMN3
## AAACCCATCCTAGCCT-1_1 -0.028925620 -0.0132582758 G1 Nb5d.INF rep2 EMN1
## AAACCCATCGAAACAA-1_1 -0.008077289 -0.0028336129 G1 Nb5d.PBS rep3 IPAN4
## AAACCCATCGGTCAGC-1_1 -0.023612751 0.0327239644 G2M Nb5d.PBS rep4 EMN1
## sample tissue nCount_SCT nFeature_SCT condition
## AAACCCACAAGACGAC-1_1 Nb5d.INF4 Ileum 2598 1771 INF_CTL
## AAACCCAGTGGGCTCT-1_1 Nb5d.PBS4 Ileum 1696 971 PBS_CTL
## AAACCCAGTTTGTTCT-1_1 Nb5d.INF1 Ileum 2493 1558 INF_CTL
## AAACCCATCCTAGCCT-1_1 Nb5d.INF2 Ileum 2319 1427 INF_CTL
## AAACCCATCGAAACAA-1_1 Nb5d.PBS3 Ileum 2543 1620 PBS_CTL
## AAACCCATCGGTCAGC-1_1 Nb5d.PBS4 Ileum 2166 1274 PBS_CTL
## integrated_snn_res.2 seurat_clusters
## AAACCCACAAGACGAC-1_1 8 8
## AAACCCAGTGGGCTCT-1_1 7 7
## AAACCCAGTTTGTTCT-1_1 8 8
## AAACCCATCCTAGCCT-1_1 2 2
## AAACCCATCGAAACAA-1_1 18 18
## AAACCCATCGGTCAGC-1_1 2 2DimPlot(GEX.seur, reduction = "umap", label = T, label.size = 3.25, group.by = "newAnno") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
DimPlot(GEX.seur, reduction = "umap", label = T, label.size = 3.75, group.by = "seurat_clusters") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", split.by = "cnt", ncol = 4, cols = color.cnt)
DimPlot(subset(GEX.seur, subset = cnt %in% c("Nb5d.PBS","Nb5d.INF")), label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", cols = color.cnt[1:2]) +
DimPlot(subset(GEX.seur, subset = cnt %in% c("Nb5d.CTL","Nb5d.CKO")), label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", cols = color.cnt[3:4])
re-clustering
DefaultAssay(GEX.seur) <- "integrated"
GEX.seur## An object of class Seurat
## 47356 features across 19649 samples within 3 assays
## Active assay: integrated (2397 features, 2397 variable features)
## 2 other assays present: RNA, SCT
## 3 dimensional reductions calculated: pca, tsne, umapElbowPlot(GEX.seur, ndims = 100)
ElbowPlot(GEX.seur, ndims = 50)
PCsct <- 1:20
GEX.seur <- FindNeighbors(GEX.seur, k.param = 18, dims = PCsct, compute.SNN = T, reduction = 'pca', verbose = T)## Computing nearest neighbor graph## Computing SNNGEX.seur <- FindClusters(GEX.seur, dims.use = PCsct, algorithm = 1, save.SNN =T, resolution =2, reduction = 'pca', verbose = T)## Warning: The following arguments are not used: dims.use, save.SNN, reduction## Suggested parameter: dims instead of dims.use## Warning: The following arguments are not used: dims.use, save.SNN, reduction## Suggested parameter: dims instead of dims.use## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 19649
## Number of edges: 584347
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8315
## Number of communities: 31
## Elapsed time: 2 secondsDimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "seurat_clusters") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters")
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "newAnno") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "newAnno")
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "newAnno")
DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", split.by = "cnt", ncol = 4, cols = color.cnt)
DefaultAssay(GEX.seur) <- "SCT"
GEX.seur## An object of class Seurat
## 47356 features across 19649 samples within 3 assays
## Active assay: SCT (20434 features, 0 variable features)
## 2 other assays present: RNA, integrated
## 3 dimensional reductions calculated: pca, tsne, umap# sort_clusters
GEX.seur$sort_clusters <- factor(as.character(GEX.seur$seurat_clusters),
levels = c(0,1,11,5,16,19, # EMN1
15,8, # EMN2
4,10, # EMN3
17, # EMN4
25, # EMN5
6,3,2,20, # IMN1
14, # IMN2
9,28, # IMN3
26, # IMN4
12, # IN1
24, # IN2
30, # IN3
21,22,13, 29, # IPAN1
7, 18, # IPAN2
27, # IPAN3
23 # IPAN4
))
# ttAnno1 (based on individual newAnno)
# hold
# ttAnno2 (sub-cluster IPAN1/2 as .1 & .2)
# holdmarkers.new.ss <- list(EMN=c("Chat","Bnc2","Tox","Ptprt",
"Gfra2","Oprk1","Adamtsl1",
"Fbxw15","Fbxw24","Chrna7",
"Satb1","Itga6","Cntnap5b",
"Pgm5","Chgb",
"Nxph1",
"Gabrb1","Glp2r","Nebl",
"Lrrc7",
"Ryr3","Eda",
"Hgf","Lama2","Efnb2",
"Tac1",
"Kctd8",
"Ptn",
"Ntrk2","Penk","Itgb8",
"Fut9","Nfatc1","Egfr","Pdia5",
"Ahr","Mgll",
"Aff3",
"Chrm3"
),
IMN=c("Nos1","Kcnab1",
"Gfra1","Etv1",
"Man1a","Airn",
"Adcy2",
"Col25a1",
"Cmah","Creb5","Vip","Pde1a",
"Ebf1","Gpc5","Mid1","Igfbp5",
"Ppara",
"Pcdh11x","Adcy8","Grp"
),
IN=c("Npas3","Synpr","St18","Gal",
"Nova1",
"Cdh10","Kcnk13",
"Neurod6","Moxd1","Sctr",
"Piezo1","Vipr2","Adamts9","Sst","Kcnn2"
),
IPAN=c("Calb2","Adcy1","Calcb","Nmu","Adgrg6","Pcdh10",
"Ngfr","Galr1","Il7","Aff2",
"Gpr149","Itgb6","Met","Itgbl1",
"Cdh6","Cdh8",
"Clstn2","Ano2","Ntrk3",
"Cpne4","Vwc2l","Cdh9",
"Car10","Dcc",
"Scgn",
"Vcan","Cck","Piezo2","Kcnh7",
"Rerg","Bmpr1b","Skap1","Ntng1",
"Tafa2","Nxph2"),
CONTAM=c("Actl6b","Phox2b"),
pDEGs=c("Grid2","Ide","Btaf1","Slc15a2","Ccser1","Hdac9","Rspo2","Grem2"))check.markers.ss <- as.vector(unlist(markers.new.ss))
check.markers.ss## [1] "Chat" "Bnc2" "Tox" "Ptprt" "Gfra2" "Oprk1"
## [7] "Adamtsl1" "Fbxw15" "Fbxw24" "Chrna7" "Satb1" "Itga6"
## [13] "Cntnap5b" "Pgm5" "Chgb" "Nxph1" "Gabrb1" "Glp2r"
## [19] "Nebl" "Lrrc7" "Ryr3" "Eda" "Hgf" "Lama2"
## [25] "Efnb2" "Tac1" "Kctd8" "Ptn" "Ntrk2" "Penk"
## [31] "Itgb8" "Fut9" "Nfatc1" "Egfr" "Pdia5" "Ahr"
## [37] "Mgll" "Aff3" "Chrm3" "Nos1" "Kcnab1" "Gfra1"
## [43] "Etv1" "Man1a" "Airn" "Adcy2" "Col25a1" "Cmah"
## [49] "Creb5" "Vip" "Pde1a" "Ebf1" "Gpc5" "Mid1"
## [55] "Igfbp5" "Ppara" "Pcdh11x" "Adcy8" "Grp" "Npas3"
## [61] "Synpr" "St18" "Gal" "Nova1" "Cdh10" "Kcnk13"
## [67] "Neurod6" "Moxd1" "Sctr" "Piezo1" "Vipr2" "Adamts9"
## [73] "Sst" "Kcnn2" "Calb2" "Adcy1" "Calcb" "Nmu"
## [79] "Adgrg6" "Pcdh10" "Ngfr" "Galr1" "Il7" "Aff2"
## [85] "Gpr149" "Itgb6" "Met" "Itgbl1" "Cdh6" "Cdh8"
## [91] "Clstn2" "Ano2" "Ntrk3" "Cpne4" "Vwc2l" "Cdh9"
## [97] "Car10" "Dcc" "Scgn" "Vcan" "Cck" "Piezo2"
## [103] "Kcnh7" "Rerg" "Bmpr1b" "Skap1" "Ntng1" "Tafa2"
## [109] "Nxph2" "Actl6b" "Phox2b" "Grid2" "Ide" "Btaf1"
## [115] "Slc15a2" "Ccser1" "Hdac9" "Rspo2" "Grem2"length(check.markers.ss)## [1] 119sum(check.markers.ss %in% rownames(GEX.seur))## [1] 119pm.All_new.a <- DotPlot(GEX.seur, features = check.markers.ss, group.by = "sort_clusters",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All - sort_clusters")
pm.All_new.a
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "sort_clusters")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "sort_clusters")
mod1
DefaultAssay(GEX.seur) <- "integrated"
GEX.seur## An object of class Seurat
## 47356 features across 19649 samples within 3 assays
## Active assay: integrated (2397 features, 2397 variable features)
## 2 other assays present: RNA, SCT
## 3 dimensional reductions calculated: pca, tsne, umap# remove C29
GEX.seur <- subset(GEX.seur, subset= seurat_clusters %in% c(0:28,30))
GEX.seur## An object of class Seurat
## 47356 features across 19460 samples within 3 assays
## Active assay: integrated (2397 features, 2397 variable features)
## 2 other assays present: RNA, SCT
## 3 dimensional reductions calculated: pca, tsne, umapDimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "newAnno")
ElbowPlot(GEX.seur, ndims = 50)
PCsct <- 1:18
GEX.seur <- FindNeighbors(GEX.seur, k.param = 12, dims = PCsct, compute.SNN = T, reduction = 'pca', verbose = T)## Computing nearest neighbor graph## Computing SNNGEX.seur <- FindClusters(GEX.seur, dims.use = PCsct, algorithm = 1, save.SNN =T, resolution =2.5, reduction = 'pca', verbose = T)## Warning: The following arguments are not used: dims.use, save.SNN, reduction## Suggested parameter: dims instead of dims.use## Warning: The following arguments are not used: dims.use, save.SNN, reduction## Suggested parameter: dims instead of dims.use## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 19460
## Number of edges: 442022
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8093
## Number of communities: 35
## Elapsed time: 2 secondsDimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "seurat_clusters") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "seurat_clusters")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "seurat_clusters")
mod2
# remove C34
GEX.seur <- subset(GEX.seur, subset= seurat_clusters %in% c(0:33))
GEX.seur## An object of class Seurat
## 47356 features across 19418 samples within 3 assays
## Active assay: integrated (2397 features, 2397 variable features)
## 2 other assays present: RNA, SCT
## 3 dimensional reductions calculated: pca, tsne, umapDimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "newAnno")
ElbowPlot(GEX.seur, ndims = 50)
PCsct <- 1:18
GEX.seur <- FindNeighbors(GEX.seur, k.param = 12, dims = PCsct, compute.SNN = T, reduction = 'pca', verbose = T)## Computing nearest neighbor graph## Computing SNNGEX.seur <- FindClusters(GEX.seur, dims.use = PCsct, algorithm = 1, save.SNN =T, resolution =2.5, reduction = 'pca', verbose = T)## Warning: The following arguments are not used: dims.use, save.SNN, reduction## Suggested parameter: dims instead of dims.use## Warning: The following arguments are not used: dims.use, save.SNN, reduction## Suggested parameter: dims instead of dims.use## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 19418
## Number of edges: 441598
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8095
## Number of communities: 34
## Elapsed time: 2 secondspp.umap1 <- DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "seurat_clusters") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters")
pp.umap1
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0.01, group.by = "seurat_clusters")
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 2, pt.size = 0, group.by = "seurat_clusters")
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "newAnno") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "newAnno")
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "seurat_clusters") +
DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "newAnno")
DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", split.by = "cnt", ncol = 4, cols = color.cnt)
intAnno
set modified annotations for integration as intAnno1/2
DefaultAssay(GEX.seur) <- "SCT"
GEX.seur## An object of class Seurat
## 47356 features across 19418 samples within 3 assays
## Active assay: SCT (20434 features, 0 variable features)
## 2 other assays present: RNA, integrated
## 3 dimensional reductions calculated: pca, tsne, umapGEX.seur@meta.data[,grep("snn",colnames(GEX.seur@meta.data))] <- NULL
head(GEX.seur@meta.data)## orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCACAAGACGAC-1_1 Nb5d.PBS_INF 3257 1801 0.36843721 0.3991403
## AAACCCAGTGGGCTCT-1_1 Nb5d.PBS_INF 1511 997 0.66181337 0.4632694
## AAACCCAGTTTGTTCT-1_1 Nb5d.PBS_INF 2855 1577 0.98073555 0.3152364
## AAACCCATCCTAGCCT-1_1 Nb5d.PBS_INF 2433 1451 0.08220304 0.3699137
## AAACCCATCGAAACAA-1_1 Nb5d.PBS_INF 3129 1656 0.12783637 0.4474273
## AAACCCATCGGTCAGC-1_1 Nb5d.PBS_INF 2201 1294 0.22716947 0.2271695
## S.Score G2M.Score Phase cnt rep newAnno
## AAACCCACAAGACGAC-1_1 0.011590405 -0.0004169865 S Nb5d.INF rep4 EMN3
## AAACCCAGTGGGCTCT-1_1 -0.024203070 0.0012414826 G2M Nb5d.PBS rep4 IPAN1
## AAACCCAGTTTGTTCT-1_1 -0.013980476 0.0039329410 G2M Nb5d.INF rep1 EMN3
## AAACCCATCCTAGCCT-1_1 -0.028925620 -0.0132582758 G1 Nb5d.INF rep2 EMN1
## AAACCCATCGAAACAA-1_1 -0.008077289 -0.0028336129 G1 Nb5d.PBS rep3 IPAN4
## AAACCCATCGGTCAGC-1_1 -0.023612751 0.0327239644 G2M Nb5d.PBS rep4 EMN1
## sample tissue nCount_SCT nFeature_SCT condition
## AAACCCACAAGACGAC-1_1 Nb5d.INF4 Ileum 2598 1771 INF_CTL
## AAACCCAGTGGGCTCT-1_1 Nb5d.PBS4 Ileum 1696 971 PBS_CTL
## AAACCCAGTTTGTTCT-1_1 Nb5d.INF1 Ileum 2493 1558 INF_CTL
## AAACCCATCCTAGCCT-1_1 Nb5d.INF2 Ileum 2319 1427 INF_CTL
## AAACCCATCGAAACAA-1_1 Nb5d.PBS3 Ileum 2543 1620 PBS_CTL
## AAACCCATCGGTCAGC-1_1 Nb5d.PBS4 Ileum 2166 1274 PBS_CTL
## seurat_clusters sort_clusters
## AAACCCACAAGACGAC-1_1 11 10
## AAACCCAGTGGGCTCT-1_1 22 21
## AAACCCAGTTTGTTCT-1_1 11 8
## AAACCCATCCTAGCCT-1_1 4 1
## AAACCCATCGAAACAA-1_1 19 23
## AAACCCATCGGTCAGC-1_1 8 1# sort_clusters
GEX.seur$sort_clusters <- factor(as.character(GEX.seur$seurat_clusters),
levels = c(0,5,4,8,16,17,14, # EMN1
10,27,11, # EMN2
2,6, # EMN3
28, # EMN4
26, # EMN5
1,9,15,3,18, # IMN1
13,24, # IMN2
20, # IMN3
29, # IMN4
12, # IN1
25, # IN2
33, # IN3
22,31, 23,30, # IPAN1
7, 21, # IPAN2
32, # IPAN3
19 # IPAN4
))
# intAnno1 (based on individual newAnno)
GEX.seur$intAnno1 <- as.character(GEX.seur$seurat_clusters)
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(0,5,4,8,16,17,14)] <- "EMN1"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(10,27,11)] <- "EMN2"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(2,6)] <- "EMN3"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(28)] <- "EMN4"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(26)] <- "EMN5"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(1,9,15,3,18)] <- "IMN1"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(13,24)] <- "IMN2"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(20)] <- "IMN3"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(29)] <- "IMN4"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(12)] <- "IN1"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(25)] <- "IN2"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(33)] <- "IN3"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(22,31,23,30)] <- "IPAN1"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(7,21)] <- "IPAN2"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(32)] <- "IPAN3"
GEX.seur$intAnno1[GEX.seur$intAnno1 %in% c(19)] <- "IPAN4"
GEX.seur$intAnno1 <- factor(GEX.seur$intAnno1,
levels = c(paste0("EMN",1:5),
paste0("IMN",1:4),
paste0("IN",1:3),
paste0("IPAN",1:4)))
# intAnno2 (sub-cluster IPAN1/2 as .1 & .2)
GEX.seur$intAnno2 <- as.character(GEX.seur$seurat_clusters)
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(0,5,4,8,16,17,14)] <- "EMN1"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(10,27,11)] <- "EMN2"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(2,6)] <- "EMN3"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(28)] <- "EMN4"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(26)] <- "EMN5"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(1,9,15,3,18)] <- "IMN1"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(13,24)] <- "IMN2"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(20)] <- "IMN3"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(29)] <- "IMN4"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(12)] <- "IN1"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(25)] <- "IN2"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(33)] <- "IN3"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(22,31)] <- "IPAN1.1"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(23,30)] <- "IPAN1.2"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(7)] <- "IPAN2.1"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(21)] <- "IPAN2.2"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(32)] <- "IPAN3"
GEX.seur$intAnno2[GEX.seur$intAnno2 %in% c(19)] <- "IPAN4"
GEX.seur$intAnno2 <- factor(GEX.seur$intAnno2,
levels = c(paste0("EMN",1:5),
paste0("IMN",1:4),
paste0("IN",1:3),
paste0("IPAN",c(1.1,1.2,2.1,2.2,3,4))))pp.umap2 <- DimPlot(GEX.seur, label = T, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "intAnno1") +
DimPlot(GEX.seur, label = T, label.size = 2.65, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "intAnno2")
pp.umap2 
pp.umap3 <- DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "cnt", split.by = "cnt", ncol = 4, cols = color.cnt)
pp.umap3
markers.new.ss <- list(EMN=c("Chat","Bnc2","Tox","Ptprt",
"Gfra2","Oprk1","Adamtsl1",
"Fbxw15","Fbxw24","Chrna7",
"Satb1","Itga6","Cntnap5b",
"Pgm5","Chgb",
"Nxph1",
"Gabrb1","Glp2r","Nebl",
"Lrrc7",
"Ryr3","Eda",
"Hgf","Lama2","Efnb2",
"Tac1",
"Kctd8",
"Ptn",
"Ntrk2","Penk","Itgb8",
"Fut9","Nfatc1","Egfr","Pdia5",
"Ahr","Mgll",
"Aff3",
"Chrm3"
),
IMN=c("Nos1","Kcnab1",
"Gfra1","Etv1",
"Man1a","Airn",
"Adcy2",
"Col25a1",
"Cmah","Creb5","Vip","Pde1a",
"Ebf1","Gpc5","Mid1","Igfbp5",
"Ppara",
"Pcdh11x","Adcy8","Grp"
),
IN=c("Npas3","Synpr","St18","Gal",
"Nova1",
"Cdh10","Kcnk13",
"Neurod6","Moxd1","Sctr",
"Piezo1","Vipr2","Adamts9","Sst","Kcnn2"
),
IPAN=c("Calb2","Adcy1","Calcb","Nmu","Adgrg6","Pcdh10",
"Ngfr","Galr1","Il7","Aff2",
"Gpr149","Itgb6","Met","Itgbl1",
"Cdh6","Cdh8",
"Clstn2","Ano2","Ntrk3",
"Cpne4","Vwc2l","Cdh9",
"Car10","Dcc",
"Scgn",
"Vcan","Cck","Piezo2","Kcnh7",
"Rerg","Bmpr1b","Skap1","Ntng1",
"Tafa2","Nxph2"),
CONTAM=c("Actl6b","Phox2b"),
pDEGs=c("Grid2","Ide","Btaf1","Slc15a2","Ccser1","Hdac9","Rspo2","Grem2"))check.markers.ss <- as.vector(unlist(markers.new.ss))
check.markers.ss## [1] "Chat" "Bnc2" "Tox" "Ptprt" "Gfra2" "Oprk1"
## [7] "Adamtsl1" "Fbxw15" "Fbxw24" "Chrna7" "Satb1" "Itga6"
## [13] "Cntnap5b" "Pgm5" "Chgb" "Nxph1" "Gabrb1" "Glp2r"
## [19] "Nebl" "Lrrc7" "Ryr3" "Eda" "Hgf" "Lama2"
## [25] "Efnb2" "Tac1" "Kctd8" "Ptn" "Ntrk2" "Penk"
## [31] "Itgb8" "Fut9" "Nfatc1" "Egfr" "Pdia5" "Ahr"
## [37] "Mgll" "Aff3" "Chrm3" "Nos1" "Kcnab1" "Gfra1"
## [43] "Etv1" "Man1a" "Airn" "Adcy2" "Col25a1" "Cmah"
## [49] "Creb5" "Vip" "Pde1a" "Ebf1" "Gpc5" "Mid1"
## [55] "Igfbp5" "Ppara" "Pcdh11x" "Adcy8" "Grp" "Npas3"
## [61] "Synpr" "St18" "Gal" "Nova1" "Cdh10" "Kcnk13"
## [67] "Neurod6" "Moxd1" "Sctr" "Piezo1" "Vipr2" "Adamts9"
## [73] "Sst" "Kcnn2" "Calb2" "Adcy1" "Calcb" "Nmu"
## [79] "Adgrg6" "Pcdh10" "Ngfr" "Galr1" "Il7" "Aff2"
## [85] "Gpr149" "Itgb6" "Met" "Itgbl1" "Cdh6" "Cdh8"
## [91] "Clstn2" "Ano2" "Ntrk3" "Cpne4" "Vwc2l" "Cdh9"
## [97] "Car10" "Dcc" "Scgn" "Vcan" "Cck" "Piezo2"
## [103] "Kcnh7" "Rerg" "Bmpr1b" "Skap1" "Ntng1" "Tafa2"
## [109] "Nxph2" "Actl6b" "Phox2b" "Grid2" "Ide" "Btaf1"
## [115] "Slc15a2" "Ccser1" "Hdac9" "Rspo2" "Grem2"length(check.markers.ss)## [1] 119sum(check.markers.ss %in% rownames(GEX.seur))## [1] 119pm.All_new.a <- DotPlot(GEX.seur, features = check.markers.ss, group.by = "sort_clusters",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All - sort_clusters")
pm.All_new.a
pm.All_new.b <- DotPlot(GEX.seur, features = check.markers.ss, group.by = "intAnno1",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All - intAnno1")
pm.All_new.b
pm.All_new.c <- DotPlot(GEX.seur, features = check.markers.ss, group.by = "intAnno2",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All - intAnno2")
pm.All_new.c
#
pm.All_new.c1 <- DotPlot(subset(GEX.seur, subset= cnt %in% c("Nb5d.PBS")), features = check.markers.ss, group.by = "intAnno2",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="Nb5d.PBS - intAnno2")
pm.All_new.c1
#
pm.All_new.c2 <- DotPlot(subset(GEX.seur, subset= cnt %in% c("Nb5d.INF")), features = check.markers.ss, group.by = "intAnno2",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="Nb5d.INF - intAnno2")
pm.All_new.c2
#
pm.All_new.c3 <- DotPlot(subset(GEX.seur, subset= cnt %in% c("Nb5d.CTL")), features = check.markers.ss, group.by = "intAnno2",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="Nb5d.CTL - intAnno2")
pm.All_new.c3
#
pm.All_new.c4 <- DotPlot(subset(GEX.seur, subset= cnt %in% c("Nb5d.CKO")), features = check.markers.ss, group.by = "intAnno2",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="Nb5d.CKO - intAnno2")
pm.All_new.c4
#saveRDS(GEX.seur, "./sn10x_WYS.sct_anno.s.rds")
#GEX.seur <- readRDS("./sn10x_WYS.sct_anno.s.rds")cell composition
pumap.test0 <- DimPlot(GEX.seur, reduction = "umap", label = T, label.size = 2.75, group.by = "intAnno2") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
pumap.test0
pumap.test0s <- DimPlot(GEX.seur, reduction = "umap", label = T, label.size = 2.75, group.by = "intAnno1") +
DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
pumap.test0s
PBS_INF
test1.seur <- subset(GEX.seur, subset= cnt %in% c("Nb5d.PBS","Nb5d.INF"))
test1.seur## An object of class Seurat
## 47356 features across 8232 samples within 3 assays
## Active assay: SCT (20434 features, 0 variable features)
## 2 other assays present: RNA, integrated
## 3 dimensional reductions calculated: pca, tsne, umapcolor.test1 <- color.cnt[1:2]pumap.test1 <- DimPlot(test1.seur, reduction = "umap", label = T, label.size = 2.75, group.by = "intAnno2") +
DimPlot(test1.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.test1)
pumap.test1
pumap.test1s <- DimPlot(test1.seur, reduction = "umap", label = T, label.size = 2.75, group.by = "intAnno1") +
DimPlot(test1.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.test1)
pumap.test1s
heatmap
cowplot::plot_grid(
pheatmap::pheatmap(table(cnt=test1.seur$sample,
clusters=test1.seur$intAnno1)[c(5:8,1:4),],
main = "Cell Count",
gaps_row = c(4),
gaps_col = c(5,9,12),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,
pheatmap::pheatmap(100*rowRatio(table(cnt=test1.seur$sample,
clusters=test1.seur$intAnno1)[c(5:8,1:4),] ),
main = "Cell Ratio",
gaps_row = c(4),
gaps_col = c(5,9,12),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)
cowplot::plot_grid(
pheatmap::pheatmap(table(cnt=test1.seur$sample,
clusters=test1.seur$intAnno2)[c(5:8,1:4),],
main = "Cell Count",
gaps_row = c(4),
gaps_col = c(5,9,12),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,
pheatmap::pheatmap(100*rowRatio(table(cnt=test1.seur$sample,
clusters=test1.seur$intAnno2)[c(5:8,1:4),] ),
main = "Cell Ratio",
gaps_row = c(4),
gaps_col = c(5,9,12),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)
barplot
stat1_intAnno1 <- test1.seur@meta.data[,c("cnt","rep","intAnno1")]
stat1_intAnno1.s <- stat1_intAnno1 %>%
group_by(cnt,rep,intAnno1) %>%
summarise(count=n()) %>%
mutate(prop= count/sum(count)) %>%
ungroup## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.pcom1.intAnno1 <- stat1_intAnno1.s %>%
ggplot(aes(x = intAnno1, y = 100*prop, color=cnt)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
theme_classic(base_size = 15) +
scale_color_manual(values = color.test1, name="") +
labs(title="Cell Composition of Nb5d.PBS_INF, intAnno1", y = "Proportion") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
geom_point(aes(x=intAnno1, y = 100*prop, color= cnt),
position = position_dodge(0.75),
shape=16,alpha=0.75,size=2.15,
stroke=0.15, show.legend=F)
pcom1.intAnno1
stat1_intAnno2 <- test1.seur@meta.data[,c("cnt","rep","intAnno2")]
stat1_intAnno2.s <- stat1_intAnno2 %>%
group_by(cnt,rep,intAnno2) %>%
summarise(count=n()) %>%
mutate(prop= count/sum(count)) %>%
ungroup## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.pcom1.intAnno2 <- stat1_intAnno2.s %>%
ggplot(aes(x = intAnno2, y = 100*prop, color=cnt)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
theme_classic(base_size = 15) +
scale_color_manual(values = color.test1, name="") +
labs(title="Cell Composition of Nb5d.PBS_INF, intAnno2", y = "Proportion") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
geom_point(aes(x=intAnno2, y = 100*prop, color= cnt),
position = position_dodge(0.75),
shape=16,alpha=0.75,size=2.15,
stroke=0.15, show.legend=F)
pcom1.intAnno2
sig.test
glm.nb - anova.LRT
Sort.N <- c("Nb5d.PBS","Nb5d.INF")
glm.nb_anova.1.intAnno1 <- list()
for(x1 in 1:2){
for(x2 in 1:2){
if(x2>x1){
stat1_intAnno1.s_N <- stat1_intAnno1.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
total.N <- stat1_intAnno1.s_N %>%
group_by(cnt, rep) %>%
summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
rownames(total.N) <- paste0(as.character(total.N$cnt),
"_",
as.character(total.N$rep))
stat1_intAnno1.s_N$total <- total.N[paste0(stat1_intAnno1.s_N$cnt,
"_",
stat1_intAnno1.s_N$rep),"total"]
glm.nb_anova.1.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
for(AA in levels(stat1_intAnno1.s_N$intAnno1)){
glm.nb_anova.1.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
aaa <- stat1_intAnno1.s_N %>% filter(intAnno1==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat1_intAnno1.s_N %>% filter(intAnno1==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat1_intAnno1.s_N %>% filter(intAnno1==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error = function(e){
NA
})
})
})
}
glm.nb_anova.1.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.1.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]])
}
}
}glm.nb_anova.1.intAnno1_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.1.intAnno1)))
rownames(glm.nb_anova.1.intAnno1_df) <- names(glm.nb_anova.1.intAnno1)
colnames(glm.nb_anova.1.intAnno1_df) <- gsub("X","C",colnames(glm.nb_anova.1.intAnno1_df))
glm.nb_anova.1.intAnno1_df## EMN1 EMN2 EMN3 EMN4 EMN5
## Nb5d.PBSvsNb5d.INF 0.02915675 0.8569496 0.2429792 0.7228933 0.9970758
## IMN1 IMN2 IMN3 IMN4 IN1 IN2
## Nb5d.PBSvsNb5d.INF 0.01061606 0.8302173 0.9288441 0.09847105 0.669224 0.2176931
## IN3 IPAN1 IPAN2 IPAN3 IPAN4
## Nb5d.PBSvsNb5d.INF 0.6229168 0.4983735 0.4164602 0.7004949 0.6806156#options (scipen = 999)round(glm.nb_anova.1.intAnno1_df,4)## EMN1 EMN2 EMN3 EMN4 EMN5 IMN1 IMN2 IMN3
## Nb5d.PBSvsNb5d.INF 0.0292 0.8569 0.243 0.7229 0.9971 0.0106 0.8302 0.9288
## IMN4 IN1 IN2 IN3 IPAN1 IPAN2 IPAN3 IPAN4
## Nb5d.PBSvsNb5d.INF 0.0985 0.6692 0.2177 0.6229 0.4984 0.4165 0.7005 0.6806Sort.N <- c("Nb5d.PBS","Nb5d.INF")
glm.nb_anova.1.intAnno2 <- list()
for(x1 in 1:2){
for(x2 in 1:2){
if(x2>x1){
stat1_intAnno2.s_N <- stat1_intAnno2.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
total.N <- stat1_intAnno2.s_N %>%
group_by(cnt, rep) %>%
summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
rownames(total.N) <- paste0(as.character(total.N$cnt),
"_",
as.character(total.N$rep))
stat1_intAnno2.s_N$total <- total.N[paste0(stat1_intAnno2.s_N$cnt,
"_",
stat1_intAnno2.s_N$rep),"total"]
glm.nb_anova.1.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
for(AA in levels(stat1_intAnno2.s_N$intAnno2)){
glm.nb_anova.1.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
aaa <- stat1_intAnno2.s_N %>% filter(intAnno2==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat1_intAnno2.s_N %>% filter(intAnno2==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat1_intAnno2.s_N %>% filter(intAnno2==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error = function(e){
NA
})
})
})
}
glm.nb_anova.1.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.1.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]])
}
}
}glm.nb_anova.1.intAnno2_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.1.intAnno2)))
rownames(glm.nb_anova.1.intAnno2_df) <- names(glm.nb_anova.1.intAnno2)
colnames(glm.nb_anova.1.intAnno2_df) <- gsub("X","C",colnames(glm.nb_anova.1.intAnno2_df))
glm.nb_anova.1.intAnno2_df## EMN1 EMN2 EMN3 EMN4 EMN5
## Nb5d.PBSvsNb5d.INF 0.02915675 0.8569496 0.2429792 0.7228933 0.9970758
## IMN1 IMN2 IMN3 IMN4 IN1 IN2
## Nb5d.PBSvsNb5d.INF 0.01061606 0.8302173 0.9288441 0.09847105 0.669224 0.2176931
## IN3 IPAN1.1 IPAN1.2 IPAN2.1 IPAN2.2 IPAN3
## Nb5d.PBSvsNb5d.INF 0.6229168 0.8709588 0.400638 0.7016721 0.3976744 0.7004949
## IPAN4
## Nb5d.PBSvsNb5d.INF 0.6806156#options (scipen = 999)round(glm.nb_anova.1.intAnno2_df,4)## EMN1 EMN2 EMN3 EMN4 EMN5 IMN1 IMN2 IMN3
## Nb5d.PBSvsNb5d.INF 0.0292 0.8569 0.243 0.7229 0.9971 0.0106 0.8302 0.9288
## IMN4 IN1 IN2 IN3 IPAN1.1 IPAN1.2 IPAN2.1 IPAN2.2
## Nb5d.PBSvsNb5d.INF 0.0985 0.6692 0.2177 0.6229 0.871 0.4006 0.7017 0.3977
## IPAN3 IPAN4
## Nb5d.PBSvsNb5d.INF 0.7005 0.6806CTL_CKO
test2.seur <- subset(GEX.seur, subset= cnt %in% c("Nb5d.CTL","Nb5d.CKO"))
test2.seur## An object of class Seurat
## 47356 features across 11186 samples within 3 assays
## Active assay: SCT (20434 features, 0 variable features)
## 2 other assays present: RNA, integrated
## 3 dimensional reductions calculated: pca, tsne, umapcolor.test2 <- color.cnt[3:4]pumap.test2 <- DimPlot(test2.seur, reduction = "umap", label = T, label.size = 2.75, group.by = "intAnno2") +
DimPlot(test2.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.test2)
pumap.test2
pumap.test2s <- DimPlot(test2.seur, reduction = "umap", label = T, label.size = 2.75, group.by = "intAnno1") +
DimPlot(test2.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.test2)
pumap.test2s
heatmap
cowplot::plot_grid(
pheatmap::pheatmap(table(cnt=test2.seur$sample,
clusters=test2.seur$intAnno1)[c(4:7,1:3),],
main = "Cell Count",
gaps_row = c(4),
gaps_col = c(5,9,12),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,
pheatmap::pheatmap(100*rowRatio(table(cnt=test2.seur$sample,
clusters=test2.seur$intAnno1)[c(4:7,1:3),] ),
main = "Cell Ratio",
gaps_row = c(4),
gaps_col = c(5,9,12),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)
cowplot::plot_grid(
pheatmap::pheatmap(table(cnt=test2.seur$sample,
clusters=test2.seur$intAnno2)[c(4:7,1:3),],
main = "Cell Count",
gaps_row = c(4),
gaps_col = c(5,9,12),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,
pheatmap::pheatmap(100*rowRatio(table(cnt=test2.seur$sample,
clusters=test2.seur$intAnno2)[c(4:7,1:3),] ),
main = "Cell Ratio",
gaps_row = c(4),
gaps_col = c(5,9,12),
cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)
barplot
stat2_intAnno1 <- test2.seur@meta.data[,c("cnt","rep","intAnno1")]
stat2_intAnno1.s <- stat2_intAnno1 %>%
group_by(cnt,rep,intAnno1) %>%
summarise(count=n()) %>%
mutate(prop= count/sum(count)) %>%
ungroup## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.pcom2.intAnno1 <- stat2_intAnno1.s %>%
ggplot(aes(x = intAnno1, y = 100*prop, color=cnt)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
theme_classic(base_size = 15) +
scale_color_manual(values = color.test2, name="") +
labs(title="Cell Composition of Nb5d.CTL_CKO, intAnno1", y = "Proportion") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
geom_point(aes(x=intAnno1, y = 100*prop, color= cnt),
position = position_dodge(0.75),
shape=16,alpha=0.75,size=2.15,
stroke=0.15, show.legend=F)
pcom2.intAnno1
stat2_intAnno2 <- test2.seur@meta.data[,c("cnt","rep","intAnno2")]
stat2_intAnno2.s <- stat2_intAnno2 %>%
group_by(cnt,rep,intAnno2) %>%
summarise(count=n()) %>%
mutate(prop= count/sum(count)) %>%
ungroup## `summarise()` has grouped output by 'cnt', 'rep'. You can override using the
## `.groups` argument.pcom2.intAnno2 <- stat2_intAnno2.s %>%
ggplot(aes(x = intAnno2, y = 100*prop, color=cnt)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
theme_classic(base_size = 15) +
scale_color_manual(values = color.test2, name="") +
labs(title="Cell Composition of Nb5d.CTL_CKO, intAnno2", y = "Proportion") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
geom_point(aes(x=intAnno2, y = 100*prop, color= cnt),
position = position_dodge(0.75),
shape=16,alpha=0.75,size=2.15,
stroke=0.15, show.legend=F)
pcom2.intAnno2
sig.test
glm.nb - anova.LRT
Sort.N <- c("Nb5d.CTL","Nb5d.CKO")
glm.nb_anova.2.intAnno1 <- list()
for(x1 in 1:2){
for(x2 in 1:2){
if(x2>x1){
stat2_intAnno1.s_N <- stat2_intAnno1.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
total.N <- stat2_intAnno1.s_N %>%
group_by(cnt, rep) %>%
summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
rownames(total.N) <- paste0(as.character(total.N$cnt),
"_",
as.character(total.N$rep))
stat2_intAnno1.s_N$total <- total.N[paste0(stat2_intAnno1.s_N$cnt,
"_",
stat2_intAnno1.s_N$rep),"total"]
glm.nb_anova.2.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
for(AA in levels(stat2_intAnno1.s_N$intAnno1)){
glm.nb_anova.2.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
aaa <- stat2_intAnno1.s_N %>% filter(intAnno1==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat2_intAnno1.s_N %>% filter(intAnno1==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat2_intAnno1.s_N %>% filter(intAnno1==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error = function(e){
NA
})
})
})
}
glm.nb_anova.2.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.2.intAnno1[[paste0(Sort.N[x1],"vs",Sort.N[x2])]])
}
}
}glm.nb_anova.2.intAnno1_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.2.intAnno1)))
rownames(glm.nb_anova.2.intAnno1_df) <- names(glm.nb_anova.2.intAnno1)
colnames(glm.nb_anova.2.intAnno1_df) <- gsub("X","C",colnames(glm.nb_anova.2.intAnno1_df))
glm.nb_anova.2.intAnno1_df## EMN1 EMN2 EMN3 EMN4 EMN5
## Nb5d.CTLvsNb5d.CKO 0.006542252 0.5432009 0.9819977 0.02373406 0.3750407
## IMN1 IMN2 IMN3 IMN4 IN1 IN2
## Nb5d.CTLvsNb5d.CKO 0.4678658 0.9067255 0.6576028 0.7693935 0.7436069 0.3062167
## IN3 IPAN1 IPAN2 IPAN3 IPAN4
## Nb5d.CTLvsNb5d.CKO 0.7938345 0.7798518 0.3699458 0.4239918 0.7177276#options (scipen = 999)round(glm.nb_anova.2.intAnno1_df,4)## EMN1 EMN2 EMN3 EMN4 EMN5 IMN1 IMN2 IMN3 IMN4
## Nb5d.CTLvsNb5d.CKO 0.0065 0.5432 0.982 0.0237 0.375 0.4679 0.9067 0.6576 0.7694
## IN1 IN2 IN3 IPAN1 IPAN2 IPAN3 IPAN4
## Nb5d.CTLvsNb5d.CKO 0.7436 0.3062 0.7938 0.7799 0.3699 0.424 0.7177Sort.N <- c("Nb5d.CTL","Nb5d.CKO")
glm.nb_anova.2.intAnno2 <- list()
for(x1 in 1:2){
for(x2 in 1:2){
if(x2>x1){
stat2_intAnno2.s_N <- stat2_intAnno2.s %>% filter(cnt %in% c(Sort.N[x1],Sort.N[x2]))
total.N <- stat2_intAnno2.s_N %>%
group_by(cnt, rep) %>%
summarise(total=sum(count)) %>% ungroup() %>% as.data.frame()
rownames(total.N) <- paste0(as.character(total.N$cnt),
"_",
as.character(total.N$rep))
stat2_intAnno2.s_N$total <- total.N[paste0(stat2_intAnno2.s_N$cnt,
"_",
stat2_intAnno2.s_N$rep),"total"]
glm.nb_anova.2.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- list()
for(AA in levels(stat2_intAnno2.s_N$intAnno2)){
glm.nb_anova.2.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]][[AA]] <- tryCatch({
aaa <- stat2_intAnno2.s_N %>% filter(intAnno2==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=1000)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat2_intAnno2.s_N %>% filter(intAnno2==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=100)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error=function(e){
tryCatch({
aaa <- stat2_intAnno2.s_N %>% filter(intAnno2==AA)
aaa.glmnb <- MASS::glm.nb(count ~ cnt + offset(log(total)), data = aaa, maxit=10)
aaa.anova <- anova(aaa.glmnb, test = "LRT")
aaa.anova$'Pr(>Chi)'[2]
}, error = function(e){
NA
})
})
})
}
glm.nb_anova.2.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]] <- unlist(glm.nb_anova.2.intAnno2[[paste0(Sort.N[x1],"vs",Sort.N[x2])]])
}
}
}glm.nb_anova.2.intAnno2_df <- t(data.frame(Reduce(function(x,y){rbind(x,y)},glm.nb_anova.2.intAnno2)))
rownames(glm.nb_anova.2.intAnno2_df) <- names(glm.nb_anova.2.intAnno2)
colnames(glm.nb_anova.2.intAnno2_df) <- gsub("X","C",colnames(glm.nb_anova.2.intAnno2_df))
glm.nb_anova.2.intAnno2_df## EMN1 EMN2 EMN3 EMN4 EMN5
## Nb5d.CTLvsNb5d.CKO 0.006542252 0.5432009 0.9819977 0.02373406 0.3750407
## IMN1 IMN2 IMN3 IMN4 IN1 IN2
## Nb5d.CTLvsNb5d.CKO 0.4678658 0.9067255 0.6576028 0.7693935 0.7436069 0.3062167
## IN3 IPAN1.1 IPAN1.2 IPAN2.1 IPAN2.2 IPAN3
## Nb5d.CTLvsNb5d.CKO 0.7938345 0.6342571 0.9505649 0.1930516 0.668667 0.4239918
## IPAN4
## Nb5d.CTLvsNb5d.CKO 0.7177276#options (scipen = 999)round(glm.nb_anova.2.intAnno2_df,4)## EMN1 EMN2 EMN3 EMN4 EMN5 IMN1 IMN2 IMN3 IMN4
## Nb5d.CTLvsNb5d.CKO 0.0065 0.5432 0.982 0.0237 0.375 0.4679 0.9067 0.6576 0.7694
## IN1 IN2 IN3 IPAN1.1 IPAN1.2 IPAN2.1 IPAN2.2 IPAN3
## Nb5d.CTLvsNb5d.CKO 0.7436 0.3062 0.7938 0.6343 0.9506 0.1931 0.6687 0.424
## IPAN4
## Nb5d.CTLvsNb5d.CKO 0.7177all markers
intAnno1
# find markers for every cluster compared to all remaining cells
Idents(GEX.seur) <- "intAnno1"
GEX.seur <- PrepSCTFindMarkers(GEX.seur, assay = "SCT")#GEX.markers.pre <- FindAllMarkers(GEX.seur, only.pos = TRUE, min.pct = 0.01,
# assay = "SCT",
# test.use = "MAST",
# logfc.threshold = 0.1)
GEX.markers.pre <- read.table("Baf53cre_Nb.markers.SCT_intAnno1.202402.csv", header = TRUE, sep = ",")
#GEX.markers.pre %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)GEX.markers.pre$cluster <- factor(as.character(GEX.markers.pre$cluster),
levels = levels(GEX.seur$intAnno1))
markers.pre_t60 <- (GEX.markers.pre %>% group_by(cluster) %>%
filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-|Hsp",gene,invert = T,value = T)) %>%
top_n(n = 60, wt = avg_log2FC) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$gene
markers.pre_t120 <- (GEX.markers.pre %>% group_by(cluster) %>%
filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-|Hsp",gene,invert = T,value = T)) %>%
top_n(n = 120, wt = avg_log2FC) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$genettt = 471
ttt/60## [1] 7.85ttt/64## [1] 7.359375ttt/65## [1] 7.246154pp.t60 <- list()
for(i in 1:8){
pp.t60[[i]] <- DotPlot(GEX.seur, features = rev(markers.pre_t60[(60*i-59):(60*i)])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}
pp.t60## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
ttt = 840
ttt/60## [1] 14ttt/64## [1] 13.125ttt/65## [1] 12.92308pp.t120 <- list()
for(i in 1:14){
pp.t120[[i]] <- DotPlot(GEX.seur, features = rev(markers.pre_t120[(60*i-59):(60*i)])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}
pp.t120## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
##
## [[9]]
##
## [[10]]
##
## [[11]]
##
## [[12]]
##
## [[13]]
##
## [[14]]
intAnno2
# find markers for every cluster compared to all remaining cells
Idents(GEX.seur) <- "intAnno2"
#GEX.seur <- PrepSCTFindMarkers(GEX.seur, assay = "SCT")#GEX.markers.pre <- FindAllMarkers(GEX.seur, only.pos = TRUE, min.pct = 0.05,
# assay = "SCT",
# test.use = "MAST",
# logfc.threshold = 0.25)
GEX.markers.pre <- read.table("Baf53cre_Nb.markers.SCT_intAnno2.202402.csv", header = TRUE, sep = ",")
#GEX.markers.pre %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)GEX.markers.pre$cluster <- factor(as.character(GEX.markers.pre$cluster),
levels = levels(GEX.seur$intAnno2))
markers.pre_t60 <- (GEX.markers.pre %>% group_by(cluster) %>%
filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-|Hsp",gene,invert = T,value = T)) %>%
top_n(n = 60, wt = avg_log2FC) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$gene
markers.pre_t120 <- (GEX.markers.pre %>% group_by(cluster) %>%
filter(pct.1>0.05 & gene %in% grep("Rps|Rpl|mt-|Hsp",gene,invert = T,value = T)) %>%
top_n(n = 120, wt = avg_log2FC) %>%
ungroup() %>%
arrange(desc(avg_log2FC*pct.1),gene) %>%
distinct(gene, .keep_all = TRUE) %>%
arrange(cluster,p_val_adj))$genettt = 485
ttt/60## [1] 8.083333ttt/64## [1] 7.578125ttt/65## [1] 7.461538pp.t60 <- list()
for(i in 1:8){
pp.t60[[i]] <- DotPlot(GEX.seur, features = rev(markers.pre_t60[(64*i-63):(64*i)])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}
pp.t60## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
ttt = 875
ttt/60## [1] 14.58333ttt/64## [1] 13.67188ttt/65## [1] 13.46154pp.t120 <- list()
for(i in 1:14){
pp.t120[[i]] <- DotPlot(GEX.seur, features = rev(markers.pre_t120[(64*i-63):(64*i)])) + coord_flip() +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))
}
pp.t120## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
##
## [[9]]
##
## [[10]]
##
## [[11]]
##
## [[12]]
##
## [[13]]
##
## [[14]]
DEGs
head(GEX.seur@meta.data)## orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCACAAGACGAC-1_1 Nb5d.PBS_INF 3257 1801 0.36843721 0.3991403
## AAACCCAGTGGGCTCT-1_1 Nb5d.PBS_INF 1511 997 0.66181337 0.4632694
## AAACCCAGTTTGTTCT-1_1 Nb5d.PBS_INF 2855 1577 0.98073555 0.3152364
## AAACCCATCCTAGCCT-1_1 Nb5d.PBS_INF 2433 1451 0.08220304 0.3699137
## AAACCCATCGAAACAA-1_1 Nb5d.PBS_INF 3129 1656 0.12783637 0.4474273
## AAACCCATCGGTCAGC-1_1 Nb5d.PBS_INF 2201 1294 0.22716947 0.2271695
## S.Score G2M.Score Phase cnt rep newAnno
## AAACCCACAAGACGAC-1_1 0.011590405 -0.0004169865 S Nb5d.INF rep4 EMN3
## AAACCCAGTGGGCTCT-1_1 -0.024203070 0.0012414826 G2M Nb5d.PBS rep4 IPAN1
## AAACCCAGTTTGTTCT-1_1 -0.013980476 0.0039329410 G2M Nb5d.INF rep1 EMN3
## AAACCCATCCTAGCCT-1_1 -0.028925620 -0.0132582758 G1 Nb5d.INF rep2 EMN1
## AAACCCATCGAAACAA-1_1 -0.008077289 -0.0028336129 G1 Nb5d.PBS rep3 IPAN4
## AAACCCATCGGTCAGC-1_1 -0.023612751 0.0327239644 G2M Nb5d.PBS rep4 EMN1
## sample tissue nCount_SCT nFeature_SCT condition
## AAACCCACAAGACGAC-1_1 Nb5d.INF4 Ileum 2592 1794 INF_CTL
## AAACCCAGTGGGCTCT-1_1 Nb5d.PBS4 Ileum 1694 996 PBS_CTL
## AAACCCAGTTTGTTCT-1_1 Nb5d.INF1 Ileum 2495 1576 INF_CTL
## AAACCCATCCTAGCCT-1_1 Nb5d.INF2 Ileum 2324 1451 INF_CTL
## AAACCCATCGAAACAA-1_1 Nb5d.PBS3 Ileum 2552 1646 PBS_CTL
## AAACCCATCGGTCAGC-1_1 Nb5d.PBS4 Ileum 2171 1293 PBS_CTL
## seurat_clusters sort_clusters intAnno1 intAnno2
## AAACCCACAAGACGAC-1_1 11 11 EMN2 EMN2
## AAACCCAGTGGGCTCT-1_1 22 22 IPAN1 IPAN1.1
## AAACCCAGTTTGTTCT-1_1 11 11 EMN2 EMN2
## AAACCCATCCTAGCCT-1_1 4 4 EMN1 EMN1
## AAACCCATCGAAACAA-1_1 19 19 IPAN4 IPAN4
## AAACCCATCGGTCAGC-1_1 8 8 EMN1 EMN1Idents(GEX.seur) <- "cnt"
#GEX.seur <- PrepSCTFindMarkers(GEX.seur, assay = "SCT")neur.clusters <- grep("EMN|IMN|IN|IPAN",levels(GEX.seur$intAnno2), value = T)
neur.clusters## [1] "EMN1" "EMN2" "EMN3" "EMN4" "EMN5" "IMN1" "IMN2"
## [8] "IMN3" "IMN4" "IN1" "IN2" "IN3" "IPAN1.1" "IPAN1.2"
## [15] "IPAN2.1" "IPAN2.2" "IPAN3" "IPAN4"#
all_new <- neur.clusters
all_new## [1] "EMN1" "EMN2" "EMN3" "EMN4" "EMN5" "IMN1" "IMN2"
## [8] "IMN3" "IMN4" "IN1" "IN2" "IN3" "IPAN1.1" "IPAN1.2"
## [15] "IPAN2.1" "IPAN2.2" "IPAN3" "IPAN4"#
list_new <- list()
list_new[["All"]] <- all_new
list_new[["EMNs"]] <- grep("EMN",all_new,value = T)
list_new[["IMNs"]] <- grep("IMN",all_new,value = T)
list_new[["IPAN1"]] <- grep("IPAN1",all_new,value = T)
list_new[["IPAN2"]] <- grep("IPAN2",all_new,value = T)
#list_new[["INs"]] <- grep("IN",all_new,value = T)
#list_new[["IPANs"]] <- grep("IPAN",all_new,value = T)
for(NN in all_new){
list_new[[NN]] <- NN
}
names_new <- names(list_new)
list_new## $All
## [1] "EMN1" "EMN2" "EMN3" "EMN4" "EMN5" "IMN1" "IMN2"
## [8] "IMN3" "IMN4" "IN1" "IN2" "IN3" "IPAN1.1" "IPAN1.2"
## [15] "IPAN2.1" "IPAN2.2" "IPAN3" "IPAN4"
##
## $EMNs
## [1] "EMN1" "EMN2" "EMN3" "EMN4" "EMN5"
##
## $IMNs
## [1] "IMN1" "IMN2" "IMN3" "IMN4"
##
## $IPAN1
## [1] "IPAN1.1" "IPAN1.2"
##
## $IPAN2
## [1] "IPAN2.1" "IPAN2.2"
##
## $EMN1
## [1] "EMN1"
##
## $EMN2
## [1] "EMN2"
##
## $EMN3
## [1] "EMN3"
##
## $EMN4
## [1] "EMN4"
##
## $EMN5
## [1] "EMN5"
##
## $IMN1
## [1] "IMN1"
##
## $IMN2
## [1] "IMN2"
##
## $IMN3
## [1] "IMN3"
##
## $IMN4
## [1] "IMN4"
##
## $IN1
## [1] "IN1"
##
## $IN2
## [1] "IN2"
##
## $IN3
## [1] "IN3"
##
## $IPAN1.1
## [1] "IPAN1.1"
##
## $IPAN1.2
## [1] "IPAN1.2"
##
## $IPAN2.1
## [1] "IPAN2.1"
##
## $IPAN2.2
## [1] "IPAN2.2"
##
## $IPAN3
## [1] "IPAN3"
##
## $IPAN4
## [1] "IPAN4"names_new## [1] "All" "EMNs" "IMNs" "IPAN1" "IPAN2" "EMN1" "EMN2"
## [8] "EMN3" "EMN4" "EMN5" "IMN1" "IMN2" "IMN3" "IMN4"
## [15] "IN1" "IN2" "IN3" "IPAN1.1" "IPAN1.2" "IPAN2.1" "IPAN2.2"
## [22] "IPAN3" "IPAN4"PBS_INF
test1.seur <- subset(GEX.seur, subset= cnt %in% c("Nb5d.PBS","Nb5d.INF"))
test1.seur## An object of class Seurat
## 47356 features across 8232 samples within 3 assays
## Active assay: SCT (20434 features, 0 variable features)
## 2 other assays present: RNA, integrated
## 3 dimensional reductions calculated: pca, tsne, umaprun
#test1.DEGs_new# save DEGs' table
df_test1.DEGs_new <- data.frame()
for(nn in names_new){
df_test1.DEGs_new <- rbind(df_test1.DEGs_new, data.frame(test1.DEGs_new[[nn]],intAnno2=nn))
}
#write.csv(df_test1.DEGs_new, "Baf53cre_Nb.DEGs.PBSvsINF.intAnno2.csv")stat
df_test1.DEGs_new$intAnno2 <- factor(as.character(df_test1.DEGs_new$intAnno2),
levels = names_new)cut.padj = 0.05
cut.log2FC = 0.1
cut.pct1 = 0.05
stat_test1a.DEGs_new <- df_test1.DEGs_new %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(intAnno2,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame()
stat_test1a.DEGs_new## intAnno2 cluster up.DEGs
## 1 All Nb5d.PBS 104
## 2 All Nb5d.INF 128
## 3 EMNs Nb5d.PBS 47
## 4 EMNs Nb5d.INF 56
## 5 IMNs Nb5d.PBS 11
## 6 IMNs Nb5d.INF 14
## 7 IPAN1 Nb5d.PBS 43
## 8 IPAN1 Nb5d.INF 127
## 9 IPAN2 Nb5d.PBS 15
## 10 IPAN2 Nb5d.INF 63
## 11 EMN1 Nb5d.PBS 21
## 12 EMN1 Nb5d.INF 31
## 13 EMN2 Nb5d.PBS 6
## 14 EMN2 Nb5d.INF 14
## 15 EMN3 Nb5d.INF 4
## 16 EMN4 Nb5d.INF 2
## 17 IMN1 Nb5d.PBS 3
## 18 IMN1 Nb5d.INF 5
## 19 IMN2 Nb5d.INF 1
## 20 IMN3 Nb5d.PBS 1
## 21 IMN3 Nb5d.INF 1
## 22 IMN4 Nb5d.PBS 1
## 23 IMN4 Nb5d.INF 1
## 24 IN1 Nb5d.INF 1
## 25 IPAN1.1 Nb5d.PBS 17
## 26 IPAN1.1 Nb5d.INF 71
## 27 IPAN1.2 Nb5d.PBS 22
## 28 IPAN1.2 Nb5d.INF 74
## 29 IPAN2.1 Nb5d.PBS 9
## 30 IPAN2.1 Nb5d.INF 60
## 31 IPAN2.2 Nb5d.INF 7
## 32 IPAN4 Nb5d.INF 5df_test1.DEGs_new %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(intAnno2,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame() %>%
ggplot(aes(x=intAnno2, y=up.DEGs, color = cluster)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
theme_classic(base_size = 15) +
scale_color_manual(values = color.test1, name="") +
labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Count") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
title =element_text(size=12, face='bold'))
cut.padj = 0.05
cut.log2FC = 0.2
cut.pct1 = 0.05
stat_test1b.DEGs_new <- df_test1.DEGs_new %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(intAnno2,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame()
stat_test1b.DEGs_new## intAnno2 cluster up.DEGs
## 1 All Nb5d.PBS 73
## 2 All Nb5d.INF 98
## 3 EMNs Nb5d.PBS 41
## 4 EMNs Nb5d.INF 53
## 5 IMNs Nb5d.PBS 10
## 6 IMNs Nb5d.INF 13
## 7 IPAN1 Nb5d.PBS 43
## 8 IPAN1 Nb5d.INF 127
## 9 IPAN2 Nb5d.PBS 15
## 10 IPAN2 Nb5d.INF 63
## 11 EMN1 Nb5d.PBS 18
## 12 EMN1 Nb5d.INF 29
## 13 EMN2 Nb5d.PBS 6
## 14 EMN2 Nb5d.INF 14
## 15 EMN3 Nb5d.INF 4
## 16 EMN4 Nb5d.INF 2
## 17 IMN1 Nb5d.PBS 3
## 18 IMN1 Nb5d.INF 5
## 19 IMN2 Nb5d.INF 1
## 20 IMN3 Nb5d.PBS 1
## 21 IMN3 Nb5d.INF 1
## 22 IMN4 Nb5d.PBS 1
## 23 IMN4 Nb5d.INF 1
## 24 IN1 Nb5d.INF 1
## 25 IPAN1.1 Nb5d.PBS 17
## 26 IPAN1.1 Nb5d.INF 71
## 27 IPAN1.2 Nb5d.PBS 22
## 28 IPAN1.2 Nb5d.INF 74
## 29 IPAN2.1 Nb5d.PBS 9
## 30 IPAN2.1 Nb5d.INF 60
## 31 IPAN2.2 Nb5d.INF 7
## 32 IPAN4 Nb5d.INF 5df_test1.DEGs_new %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(intAnno2,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame() %>%
ggplot(aes(x=intAnno2, y=up.DEGs, color = cluster)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
theme_classic(base_size = 15) +
scale_color_manual(values = color.test1, name="") +
labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Count") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
title =element_text(size=12, face='bold'))
plot
#save(test1.DEGs_new, file = "test1.DEGs_new.Rdata")pp_test1.DEGs.new <- lapply(list_new, function(x){NA})#
test1.DEGs_new.combine <- test1.DEGs_new
test1.DEGs_new.combine <- lapply(test1.DEGs_new.combine, function(x){
x[x$cluster=="Nb5d.PBS","avg_log2FC"] <- -x[x$cluster=="Nb5d.PBS","avg_log2FC"]
x
})All
pp_test1.DEGs.new$All <- test1.DEGs_new.combine$All %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="All Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$All
EMNs
pp_test1.DEGs.new$EMNs <- test1.DEGs_new.combine$EMNs %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="EMNs Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$EMNs
EMN1
pp_test1.DEGs.new$EMN1 <- test1.DEGs_new.combine$EMN1 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="EMN1 Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$EMN1
EMN2
pp_test1.DEGs.new$EMN2 <- test1.DEGs_new.combine$EMN2 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="EMN2 Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$EMN2
EMN3
pp_test1.DEGs.new$EMN3 <- test1.DEGs_new.combine$EMN3 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="EMN3 Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$EMN3
EMN4
pp_test1.DEGs.new$EMN4 <- test1.DEGs_new.combine$EMN4 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="EMN4 Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$EMN4
IMNs
pp_test1.DEGs.new$IMNs <- test1.DEGs_new.combine$IMNs %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="IMNs Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IMNs
IMN1
pp_test1.DEGs.new$IMN1 <- test1.DEGs_new.combine$IMN1 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="IMN1 Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IMN1
IMN2
pp_test1.DEGs.new$IMN2 <- test1.DEGs_new.combine$IMN2 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="IMN2 Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IMN2
IMN3
pp_test1.DEGs.new$IMN3 <- test1.DEGs_new.combine$IMN3 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="IMN3 Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IMN3
IMN4
pp_test1.DEGs.new$IMN4 <- test1.DEGs_new.combine$IMN4 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="IMN4 Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IMN4
IN1
pp_test1.DEGs.new$IN1 <- test1.DEGs_new.combine$IN1 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="IN1 Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IN1
IPAN1
pp_test1.DEGs.new$IPAN1 <- test1.DEGs_new.combine$IPAN1 %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="IPAN1 Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN1
IPAN1.1
pp_test1.DEGs.new$IPAN1.1 <- test1.DEGs_new.combine$IPAN1.1 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="IPAN1.1 Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN1.1
IPAN1.2
pp_test1.DEGs.new$IPAN1.2 <- test1.DEGs_new.combine$IPAN1.2 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="IPAN1.2 Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN1.2
IPAN2
pp_test1.DEGs.new$IPAN2 <- test1.DEGs_new.combine$IPAN2 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="IPAN2 Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN2
IPAN2.1
pp_test1.DEGs.new$IPAN2.1 <- test1.DEGs_new.combine$IPAN2.1 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="IPAN2.1 Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN2.1
IPAN2.2
pp_test1.DEGs.new$IPAN2.2 <- test1.DEGs_new.combine$IPAN2.2 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="IPAN2.2 Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN2.2
IPAN4
pp_test1.DEGs.new$IPAN4 <- test1.DEGs_new.combine$IPAN4 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.INF","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.PBS"=color.test1[1],
"Nb5d.INF"=color.test1[2],
"None"="grey")) +
theme_classic() + labs(title="IPAN4 Nb5d.PBS vs Nb5d.INF") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test1.DEGs.new$IPAN4
neur.ol <- c("All","EMNs","EMN1","EMN2","IMNs","IMN1","IPAN1","IPAN1.1","IPAN1.2","IPAN2","IPAN2.1")
neur.ol## [1] "All" "EMNs" "EMN1" "EMN2" "IMNs" "IMN1" "IPAN1"
## [8] "IPAN1.1" "IPAN1.2" "IPAN2" "IPAN2.1"CTL_CKO
test2.seur <- subset(GEX.seur, subset= cnt %in% c("Nb5d.CTL","Nb5d.CKO"))
test2.seur## An object of class Seurat
## 47356 features across 11186 samples within 3 assays
## Active assay: SCT (20434 features, 0 variable features)
## 2 other assays present: RNA, integrated
## 3 dimensional reductions calculated: pca, tsne, umaprun
#test2.DEGs_new# save DEGs' table
df_test2.DEGs_new <- data.frame()
for(nn in names_new){
df_test2.DEGs_new <- rbind(df_test2.DEGs_new, data.frame(test2.DEGs_new[[nn]],intAnno2=nn))
}
#write.csv(df_test2.DEGs_new, "Baf53cre_Nb.DEGs.CTLvsCKO.intAnno2.csv")stat
df_test2.DEGs_new$intAnno2 <- factor(as.character(df_test2.DEGs_new$intAnno2),
levels = names_new)cut.padj = 0.05
cut.log2FC = 0.1
cut.pct1 = 0.05
stat_test2a.DEGs_new <- df_test2.DEGs_new %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(intAnno2,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame()
stat_test2a.DEGs_new## intAnno2 cluster up.DEGs
## 1 All Nb5d.CTL 159
## 2 All Nb5d.CKO 52
## 3 EMNs Nb5d.CTL 66
## 4 EMNs Nb5d.CKO 20
## 5 IMNs Nb5d.CTL 13
## 6 IMNs Nb5d.CKO 4
## 7 IPAN1 Nb5d.CTL 107
## 8 IPAN1 Nb5d.CKO 30
## 9 IPAN2 Nb5d.CTL 39
## 10 IPAN2 Nb5d.CKO 7
## 11 EMN1 Nb5d.CTL 22
## 12 EMN1 Nb5d.CKO 3
## 13 EMN2 Nb5d.CTL 4
## 14 EMN2 Nb5d.CKO 2
## 15 EMN3 Nb5d.CTL 5
## 16 IMN1 Nb5d.CTL 9
## 17 IMN1 Nb5d.CKO 1
## 18 IN1 Nb5d.CTL 1
## 19 IPAN1.1 Nb5d.CTL 46
## 20 IPAN1.1 Nb5d.CKO 17
## 21 IPAN1.2 Nb5d.CTL 60
## 22 IPAN1.2 Nb5d.CKO 16
## 23 IPAN2.1 Nb5d.CTL 38
## 24 IPAN2.1 Nb5d.CKO 9
## 25 IPAN2.2 Nb5d.CTL 1df_test2.DEGs_new %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(intAnno2,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame() %>%
ggplot(aes(x=intAnno2, y=up.DEGs, color = cluster)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
theme_classic(base_size = 15) +
scale_color_manual(values = color.test2, name="") +
labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Count") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
title =element_text(size=12, face='bold'))
cut.padj = 0.05
cut.log2FC = 0.2
cut.pct1 = 0.05
stat_test2b.DEGs_new <- df_test2.DEGs_new %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(intAnno2,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame()
stat_test2b.DEGs_new## intAnno2 cluster up.DEGs
## 1 All Nb5d.CTL 74
## 2 All Nb5d.CKO 8
## 3 EMNs Nb5d.CTL 59
## 4 EMNs Nb5d.CKO 9
## 5 IMNs Nb5d.CTL 13
## 6 IMNs Nb5d.CKO 1
## 7 IPAN1 Nb5d.CTL 107
## 8 IPAN1 Nb5d.CKO 30
## 9 IPAN2 Nb5d.CTL 39
## 10 IPAN2 Nb5d.CKO 7
## 11 EMN1 Nb5d.CTL 21
## 12 EMN1 Nb5d.CKO 2
## 13 EMN2 Nb5d.CTL 4
## 14 EMN2 Nb5d.CKO 1
## 15 EMN3 Nb5d.CTL 5
## 16 IMN1 Nb5d.CTL 9
## 17 IMN1 Nb5d.CKO 1
## 18 IN1 Nb5d.CTL 1
## 19 IPAN1.1 Nb5d.CTL 46
## 20 IPAN1.1 Nb5d.CKO 17
## 21 IPAN1.2 Nb5d.CTL 60
## 22 IPAN1.2 Nb5d.CKO 16
## 23 IPAN2.1 Nb5d.CTL 38
## 24 IPAN2.1 Nb5d.CKO 9
## 25 IPAN2.2 Nb5d.CTL 1df_test2.DEGs_new %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(intAnno2,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame() %>%
ggplot(aes(x=intAnno2, y=up.DEGs, color = cluster)) +
geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
theme_classic(base_size = 15) +
scale_color_manual(values = color.test2, name="") +
labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Count") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
title =element_text(size=12, face='bold'))
plot
#save(test2.DEGs_new, file = "test2.DEGs_new.Rdata")pp_test2.DEGs.new <- lapply(list_new, function(x){NA})#
test2.DEGs_new.combine <- test2.DEGs_new
test2.DEGs_new.combine <- lapply(test2.DEGs_new.combine, function(x){
x[x$cluster=="Nb5d.CTL","avg_log2FC"] <- -x[x$cluster=="Nb5d.CTL","avg_log2FC"]
x
})All
pp_test2.DEGs.new$All <- test2.DEGs_new.combine$All %>%
mutate(label=ifelse(((p_val_adj < 1e-5 & avg_log2FC<0) | (p_val_adj < 1e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
"Nb5d.CKO"=color.test2[2],
"None"="grey")) +
theme_classic() + labs(title="All Nb5d.CTL vs Nb5d.CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$All
EMNs
pp_test2.DEGs.new$EMNs <- test2.DEGs_new.combine$EMNs %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
"Nb5d.CKO"=color.test2[2],
"None"="grey")) +
theme_classic() + labs(title="EMNs Nb5d.CTL vs Nb5d.CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$EMNs
EMN1
pp_test2.DEGs.new$EMN1 <- test2.DEGs_new.combine$EMN1 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
"Nb5d.CKO"=color.test2[2],
"None"="grey")) +
theme_classic() + labs(title="EMN1 Nb5d.CTL vs Nb5d.CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$EMN1
EMN2
pp_test2.DEGs.new$EMN2 <- test2.DEGs_new.combine$EMN2 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
"Nb5d.CKO"=color.test2[2],
"None"="grey")) +
theme_classic() + labs(title="EMN2 Nb5d.CTL vs Nb5d.CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$EMN2
EMN3
pp_test2.DEGs.new$EMN3 <- test2.DEGs_new.combine$EMN3 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
"Nb5d.CKO"=color.test2[2],
"None"="grey")) +
theme_classic() + labs(title="EMN3 Nb5d.CTL vs Nb5d.CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$EMN3
IMNs
pp_test2.DEGs.new$IMNs <- test2.DEGs_new.combine$IMNs %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
"Nb5d.CKO"=color.test2[2],
"None"="grey")) +
theme_classic() + labs(title="IMNs Nb5d.CTL vs Nb5d.CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IMNs
IMN1
pp_test2.DEGs.new$IMN1 <- test2.DEGs_new.combine$IMN1 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
"Nb5d.CKO"=color.test2[2],
"None"="grey")) +
theme_classic() + labs(title="IMN1 Nb5d.CTL vs Nb5d.CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IMN1
IN1
pp_test2.DEGs.new$IN1 <- test2.DEGs_new.combine$IN1 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
"Nb5d.CKO"=color.test2[2],
"None"="grey")) +
theme_classic() + labs(title="IN1 Nb5d.CTL vs Nb5d.CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IN1
IPAN1
pp_test2.DEGs.new$IPAN1 <- test2.DEGs_new.combine$IPAN1 %>%
mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
"Nb5d.CKO"=color.test2[2],
"None"="grey")) +
theme_classic() + labs(title="IPAN1 Nb5d.CTL vs Nb5d.CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN1
IPAN1.1
pp_test2.DEGs.new$IPAN1.1 <- test2.DEGs_new.combine$IPAN1.1 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
"Nb5d.CKO"=color.test2[2],
"None"="grey")) +
theme_classic() + labs(title="IPAN1.1 Nb5d.CTL vs Nb5d.CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN1.1
IPAN1.2
pp_test2.DEGs.new$IPAN1.2 <- test2.DEGs_new.combine$IPAN1.2 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
"Nb5d.CKO"=color.test2[2],
"None"="grey")) +
theme_classic() + labs(title="IPAN1.2 Nb5d.CTL vs Nb5d.CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN1.2
IPAN2
pp_test2.DEGs.new$IPAN2 <- test2.DEGs_new.combine$IPAN2 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
"Nb5d.CKO"=color.test2[2],
"None"="grey")) +
theme_classic() + labs(title="IPAN2 Nb5d.CTL vs Nb5d.CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN2
IPAN2.1
pp_test2.DEGs.new$IPAN2.1 <- test2.DEGs_new.combine$IPAN2.1 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
"Nb5d.CKO"=color.test2[2],
"None"="grey")) +
theme_classic() + labs(title="IPAN2.1 Nb5d.CTL vs Nb5d.CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN2.1
IPAN2.2
pp_test2.DEGs.new$IPAN2.2 <- test2.DEGs_new.combine$IPAN2.2 %>%
mutate(label=ifelse(((p_val_adj < 5e-2 & avg_log2FC<0) | (p_val_adj < 5e-2 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
#!grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
!grepl("^Rps|^Rpl",gene),gene,""),
sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"Nb5d.CTL",ifelse(p_val_adj<0.05 & avg_log2FC>0,"Nb5d.CKO","None")),
padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
geom_point(shape=21, alpha=0.65) +
geom_text_repel(size=2.5, max.overlaps = 500) +
scale_fill_manual(values = c("Nb5d.CTL"=color.test2[1],
"Nb5d.CKO"=color.test2[2],
"None"="grey")) +
theme_classic() + labs(title="IPAN2.2 Nb5d.CTL vs Nb5d.CKO") +
guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")
pp_test2.DEGs.new$IPAN2.2
signature scores
top list
use top120-built top markers
load("I:/Shared_win/projects/20231220_10x_SZJ/analysis_plus/top.markerlist.RData")data.frame(row.names = "topmarkers",lapply(top.list,length))## EMN1 EMN2 EMN3 EMN4 EMN5 IMN1 IMN2 IMN3 IMN4 IN1 IN2 IN3 IPAN1.1
## topmarkers 96 32 37 34 70 71 33 49 53 60 56 63 70
## IPAN1.2 IPAN2.1 IPAN2.2 IPAN3 IPAN4
## topmarkers 51 43 52 52 50neur.clusters## [1] "EMN1" "EMN2" "EMN3" "EMN4" "EMN5" "IMN1" "IMN2"
## [8] "IMN3" "IMN4" "IN1" "IN2" "IN3" "IPAN1.1" "IPAN1.2"
## [15] "IPAN2.1" "IPAN2.2" "IPAN3" "IPAN4"for(nn in neur.clusters){
GEX.seur <- add_geneset_score(obj = GEX.seur,
Assay = "SCT",
geneset = top.list[[nn]],
setname = paste0("score.",nn))
}## Summarizing data
## Summarizing data
## Summarizing data
## Summarizing data
## Summarizing data
## Summarizing data
## Summarizing data
## Summarizing data
## Summarizing data
## Summarizing data
## Summarizing data
## Summarizing data
## Summarizing data
## Summarizing data
## Summarizing data
## Summarizing data
## Summarizing data
## Summarizing datamapal <- colorRampPalette(RColorBrewer::brewer.pal(4,"Spectral"))(120)
FeaturePlot(GEX.seur, features = paste0("score.",neur.clusters), ncol = 5) &
scale_color_gradientn(colors = rev(mapal))
mapal <- colorRampPalette(RColorBrewer::brewer.pal(4,"Spectral"))(120)
FeaturePlot(GEX.seur, features = paste0("score.",neur.clusters), ncol = 5, raster = T, pt.size = 3.5) &
scale_color_gradientn(colors = rev(mapal))
vln.list <- list()
vln.list <- lapply(neur.clusters, function(x){
VlnPlot(GEX.seur, features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.25,0.65)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.4, 0.4, 0.45),
size=2.75
)
})
names(vln.list) <- neur.clusterscowplot::plot_grid(
plotlist = vln.list ,
ncol = 9)
overlapped DEGs
process
cut.padj = 0.05
cut.log2FC = 0.1
cut.pct1 = 0.05
stat_test1a.DEGs_new <- df_test1.DEGs_new %>% filter(p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1) %>%
group_by(intAnno2,cluster) %>%
summarise(up.DEGs = n()) %>% as.data.frame()
#stat_test1a.DEGs_newhead(df_test1.DEGs_new)## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster
## Ctnna3 8.684419e-51 0.4588045 0.846 0.761 2.129854e-46 Nb5d.PBS
## Malat1 1.335552e-36 0.1410444 1.000 1.000 3.275441e-32 Nb5d.PBS
## AY036118 7.553267e-24 0.2576808 0.788 0.750 1.852439e-19 Nb5d.PBS
## Fgfr2 2.040025e-21 0.2291577 0.838 0.769 5.003161e-17 Nb5d.PBS
## 4930447N08Rik 2.483724e-18 0.2687259 0.726 0.659 6.091333e-14 Nb5d.PBS
## Prkn 2.213112e-17 0.2764906 0.516 0.421 5.427657e-13 Nb5d.PBS
## gene intAnno2
## Ctnna3 Ctnna3 All
## Malat1 Malat1 All
## AY036118 AY036118 All
## Fgfr2 Fgfr2 All
## 4930447N08Rik 4930447N08Rik All
## Prkn Prkn All#
cut.padj = 0.05
cut.log2FC = 0.1
cut.pct1 = 0.05
#
test1_up.list <- lapply(neur.ol,function(x){NA})
names(test1_up.list) <- neur.ol
for(NN in neur.ol){
test1_up.list[[NN]] <- list()
test1_up.list[[NN]][["Nb5d.PBS"]] <- (df_test1.DEGs_new %>% filter(intAnno2 == NN &
cluster == "Nb5d.PBS" &
p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1))$gene
test1_up.list[[NN]][["Nb5d.INF"]] <- (df_test1.DEGs_new %>% filter(intAnno2 == NN &
cluster == "Nb5d.INF" &
p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1))$gene
}
#test1_up.list
#
test2_up.list <- lapply(neur.ol,function(x){NA})
names(test2_up.list) <- neur.ol
for(NN in neur.ol){
test2_up.list[[NN]] <- list()
test2_up.list[[NN]][["Nb5d.CTL"]] <- (df_test2.DEGs_new %>% filter(intAnno2 == NN &
cluster == "Nb5d.CTL" &
p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1))$gene
test2_up.list[[NN]][["Nb5d.CKO"]] <- (df_test2.DEGs_new %>% filter(intAnno2 == NN &
cluster == "Nb5d.CKO" &
p_val_adj < cut.padj &
abs(avg_log2FC) > cut.log2FC &
pct.1 > cut.pct1))$gene
}
#test2_up.listlibrary(UpSetR)#
xx=neur.ol[1]
upset(fromList(c(test1_up.list[[xx]],
test2_up.list[[xx]])),
order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)
grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))
#
xx=neur.ol[2]
upset(fromList(c(test1_up.list[[xx]],
test2_up.list[[xx]])),
order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)
grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))
#
xx=neur.ol[3]
upset(fromList(c(test1_up.list[[xx]],
test2_up.list[[xx]])),
order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)
grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))
#
xx=neur.ol[4]
upset(fromList(c(test1_up.list[[xx]],
test2_up.list[[xx]])),
order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)
grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))
#
xx=neur.ol[5]
upset(fromList(c(test1_up.list[[xx]],
test2_up.list[[xx]])),
order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)
grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))
#
xx=neur.ol[6]
upset(fromList(c(test1_up.list[[xx]],
test2_up.list[[xx]])),
order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)
grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))
#
xx=neur.ol[7]
upset(fromList(c(test1_up.list[[xx]],
test2_up.list[[xx]])),
order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)
grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))
#
xx=neur.ol[8]
upset(fromList(c(test1_up.list[[xx]],
test2_up.list[[xx]])),
order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)
grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))
#
xx=neur.ol[9]
upset(fromList(c(test1_up.list[[xx]],
test2_up.list[[xx]])),
order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)
grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))
#
xx=neur.ol[10]
upset(fromList(c(test1_up.list[[xx]],
test2_up.list[[xx]])),
order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)
grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))
#
xx=neur.ol[11]
upset(fromList(c(test1_up.list[[xx]],
test2_up.list[[xx]])),
order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)
grid::grid.text(xx,x = 0.65, y=0.95, gp=grid::gpar(fontsize=20))
save
test3_up.list <- list()
test3_up.list[["INFxCTL"]] <- list()
test3_up.list[["PBSxCKO"]] <- list()
for(NN in c("All","IPAN1","IPAN1.1","IPAN1.2","IPAN2","IPAN2.1")){
test3_up.list[["INFxCTL"]][[NN]] <- intersect(test1_up.list[[NN]][["Nb5d.INF"]],test2_up.list[[NN]][["Nb5d.CTL"]])
test3_up.list[["PBSxCKO"]][[NN]] <- intersect(test1_up.list[[NN]][["Nb5d.PBS"]],test2_up.list[[NN]][["Nb5d.CKO"]])
}
test3_up.list## $INFxCTL
## $INFxCTL$All
## [1] "Fam155a" "4930432L08Rik" "Nmu" "Parm1"
## [5] "Cacnb4" "Pcsk2" "Calcb" "Rab27b"
## [9] "Nell1" "Gm30613" "Tmcc3" "Gm21847"
## [13] "Gng2" "Hnrnpll" "Cyyr1" "Grp"
##
## $INFxCTL$IPAN1
## [1] "Parm1" "Nmu" "Cadm2" "Sgip1"
## [5] "App" "Fam155a" "Dgki" "4930432L08Rik"
## [9] "Kcnh8" "Gng2" "Cacnb4" "Hnrnpll"
## [13] "Col24a1" "Gm26871" "Alk" "Calcb"
## [17] "Asxl3" "Nell1" "Pcsk2" "Abi3bp"
## [21] "Igf1r" "Lingo2" "Dpyd" "Tll1"
## [25] "Negr1" "Dysf" "Syt9" "Gm30094"
## [29] "Scn3a" "Epha7" "Itih5" "Syt1"
## [33] "Dgkg" "Ppm1h" "Ctnnd2" "Galnt13"
## [37] "Gm12709" "1700111E14Rik" "Epha6" "Rcan3"
## [41] "Ptchd4" "Rab27b" "Gm38405" "Ppp1r12b"
## [45] "Kif5a" "Pak7" "Dner" "Tmcc3"
## [49] "Fhl1" "Dclk1" "Grp" "Gpc6"
## [53] "Galm" "Map3k3" "Gm21847" "Rock1"
## [57] "Kcna4" "Gm16158" "Slco3a1" "Gm49226"
## [61] "Rad51b" "Gucy1a2" "Slco4c1" "Klf6"
## [65] "Nell1os" "Cpne7" "Gm43391" "Pcsk2os2"
## [69] "Gm10791" "Adk" "Gm48283" "Tmem8b"
## [73] "Unc13b" "Gm21798" "Syt14" "Mapk8"
## [77] "Npr2"
##
## $INFxCTL$IPAN1.1
## [1] "Nmu" "Parm1" "Cadm2" "Fam155a"
## [5] "Sgip1" "App" "Kcnh8" "Dgki"
## [9] "4930432L08Rik" "Gng2" "Hnrnpll" "Asxl3"
## [13] "Dysf" "Dgkg" "Abi3bp" "Pcsk2"
## [17] "Cacnb4" "Calcb" "Epha6" "Lingo2"
## [21] "Epha7" "Gm30094" "Igf1r" "Ppp1r12b"
## [25] "Scn3a" "Ppm1h" "Nell1" "Rcan3"
## [29] "Fhl1" "Tmcc3" "Gm12709" "1700111E14Rik"
## [33] "Grp" "Rock1" "Cpne7" "Gucy1a2"
## [37] "Klf6"
##
## $INFxCTL$IPAN1.2
## [1] "Parm1" "Nmu" "Col24a1" "App"
## [5] "4930432L08Rik" "Alk" "Sgip1" "Cadm2"
## [9] "Fam155a" "Dgki" "Cacnb4" "Gm26871"
## [13] "Gng2" "Nell1" "Lingo2" "Hnrnpll"
## [17] "Negr1" "Tll1" "Calcb" "Rab27b"
## [21] "Kcnh8" "Gm38405" "Igf1r" "Syt9"
## [25] "Galnt13" "Gpc6" "Dpyd" "Asxl3"
## [29] "Pcsk2" "Itih5" "Scn3a" "Kif5a"
## [33] "Dclk1" "Abi3bp" "Ppm1h" "Pak7"
## [37] "Gm30094" "Gm12709" "1700111E14Rik" "Ctnnd2"
## [41] "Rad51b" "Epha7" "Ptchd4" "Kcnt2"
## [45] "Rcan3" "Tmcc3" "Epha6"
##
## $INFxCTL$IPAN2
## [1] "Fam155a" "Large1" "Gm30613" "Nrxn3" "Pcsk2" "Gm21847" "Ptprz1"
## [8] "Kcnd2" "Alk" "Calcb" "Cacnb4" "Ppm1h" "Tmcc3" "Gm21798"
## [15] "Trpc3" "Arid5b" "Nell1" "Oxr1" "Galnt13" "Gm30094" "Adcy8"
## [22] "Kcnq5" "Rgs7" "Necab1" "Dner" "Dlc1" "Hnrnpll" "Dmd"
## [29] "Antxr2" "Pkp4"
##
## $INFxCTL$IPAN2.1
## [1] "Fam155a" "Large1" "Gm30613" "Nrxn3" "Pcsk2" "Gm21847" "Cacnb4"
## [8] "Calcb" "Gm21798" "Rab27b" "Alk" "Ppm1h" "Ptprz1" "Nell1"
## [15] "Tmcc3" "Arid5b" "Oxr1" "Trpc3" "Gm30094" "Adcy8" "Kcnq5"
## [22] "Dclk1" "Galnt13" "Rgs7" "Vwc2l" "Dlc1" "Necab1" "Hnrnpll"
## [29] "Antxr2" "Dner" "Setbp1"
##
##
## $PBSxCKO
## $PBSxCKO$All
## [1] "Cdh6" "Agbl4" "Pdzrn4" "Rsrp1" "Zfp804a"
##
## $PBSxCKO$IPAN1
## [1] "Zfp804a" "Fgf13" "Cdh6" "Tafa1" "Luzp2" "Ppp3ca" "Efr3a"
## [8] "Lsamp" "Arhgap6" "Ano2" "Rab3c" "Rbfox1" "Pdzrn4" "Cntn5"
## [15] "Kcnb2" "Nme6"
##
## $PBSxCKO$IPAN1.1
## [1] "Zfp804a" "Fgf13" "Tafa1" "Lsamp" "Arhgap6"
##
## $PBSxCKO$IPAN1.2
## [1] "Cntn5" "Zfp804a" "Kcnb2" "Fgf13" "Ppp3ca" "Luzp2" "Cdh6"
## [8] "Tafa1" "Pdzrn4"
##
## $PBSxCKO$IPAN2
## [1] "Dgkb"
##
## $PBSxCKO$IPAN2.1
## [1] "Dgkb" "Cdh6"plot
length(test3_up.list$INFxCTL$IPAN1)## [1] 77vln.DEGs_INFxCTL <- list()# violin for all DEGs in INFxCTL_IPAN1
# IPAN1 only
vln.DEGs_INFxCTL[["IPAN1"]] <- lapply(test3_up.list$INFxCTL$IPAN1,function(xx){
VlnPlot(subset(GEX.seur,subset=intAnno1 %in% c("IPAN1")), features = xx, group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,4)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(1, 1, 1.5),
size=2.75
)
})
names(vln.DEGs_INFxCTL[["IPAN1"]]) <- test3_up.list$INFxCTL$IPAN1#vln.DEGs_INFxCTL[["IPAN1"]]cowplot::plot_grid(
plotlist = vln.DEGs_INFxCTL[["IPAN1"]],
ncol = 4
)
length(test3_up.list$INFxCTL$IPAN2)## [1] 30# violin for all DEGs in INFxCTL_IPAN2
# IPAN2 only
vln.DEGs_INFxCTL[["IPAN2"]] <- lapply(test3_up.list$INFxCTL$IPAN2,function(xx){
VlnPlot(subset(GEX.seur,subset=intAnno1 %in% c("IPAN2")), features = xx, group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,4)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(1, 1, 1.5),
size=2.75
)
})
names(vln.DEGs_INFxCTL[["IPAN2"]]) <- test3_up.list$INFxCTL$IPAN2
# mod one without pvalue
vln.DEGs_INFxCTL[["IPAN2"]][["Fam155a"]] <- VlnPlot(subset(GEX.seur,subset=intAnno1 %in% c("IPAN2")), features = "Fam155a", group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,4)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(3.2, 3.2, 3.8),
size=2.75
)#vln.DEGs_INFxCTL[["IPAN2"]]cowplot::plot_grid(
plotlist = vln.DEGs_INFxCTL[["IPAN2"]],
ncol = 4
)
sig.score
# scoring INFxCTL in IPAN1 and IPAN2
GEX.seur <- add_geneset_score(obj = GEX.seur,
Assay = "SCT",
geneset = test3_up.list$INFxCTL$IPAN1,
setname = "score.INFxCTL_IPAN1")## Summarizing dataGEX.seur <- add_geneset_score(obj = GEX.seur,
Assay = "SCT",
geneset = test3_up.list$INFxCTL$IPAN2,
setname = "score.INFxCTL_IPAN2")## Summarizing datavln.INFxCTL <- list()
vln.INFxCTL <- lapply(c("INFxCTL_IPAN1","INFxCTL_IPAN2"), function(x){
VlnPlot(GEX.seur, features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.25,0.65)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.4, 0.4, 0.45),
size=2.75
)
})
names(vln.INFxCTL) <- c("INFxCTL_IPAN1","INFxCTL_IPAN2")# All
cowplot::plot_grid(
plotlist = vln.INFxCTL ,
ncol = 2)
# IPAN1 only
x = "INFxCTL_IPAN1"
vln.INFxCTL_IPAN1 <- VlnPlot(subset(GEX.seur, subset=intAnno1 %in% c("IPAN1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.15,0.75)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.6, 0.6, 0.65),
size=2.75
)
vln.INFxCTL_IPAN1
# IPAN1 only
x = "INFxCTL_IPAN2"
vln.INFxCTL_IPAN2 <- VlnPlot(subset(GEX.seur, subset=intAnno1 %in% c("IPAN2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.2,1.05)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.85, 0.85, 0.95),
size=2.75
)
vln.INFxCTL_IPAN2
# scoring DEGs.All
GEX.seur <- add_geneset_score(obj = GEX.seur,
Assay = "SCT",
geneset = test1_up.list$All$Nb5d.PBS,
setname = "score.All_PBSup")## Summarizing dataGEX.seur <- add_geneset_score(obj = GEX.seur,
Assay = "SCT",
geneset = test1_up.list$All$Nb5d.INF,
setname = "score.All_INFup")## Summarizing dataGEX.seur <- add_geneset_score(obj = GEX.seur,
Assay = "SCT",
geneset = test2_up.list$All$Nb5d.CTL,
setname = "score.All_CTLup")## Summarizing dataGEX.seur <- add_geneset_score(obj = GEX.seur,
Assay = "SCT",
geneset = test2_up.list$All$Nb5d.CKO,
setname = "score.All_CKOup")## Summarizing data# scoring DEGs.IPAN1
GEX.seur <- add_geneset_score(obj = GEX.seur,
Assay = "SCT",
geneset = test1_up.list$IPAN1$Nb5d.PBS,
setname = "score.IPAN1_PBSup")## Summarizing dataGEX.seur <- add_geneset_score(obj = GEX.seur,
Assay = "SCT",
geneset = test1_up.list$IPAN1$Nb5d.INF,
setname = "score.IPAN1_INFup")## Summarizing dataGEX.seur <- add_geneset_score(obj = GEX.seur,
Assay = "SCT",
geneset = test2_up.list$IPAN1$Nb5d.CTL,
setname = "score.IPAN1_CTLup")## Summarizing dataGEX.seur <- add_geneset_score(obj = GEX.seur,
Assay = "SCT",
geneset = test2_up.list$IPAN1$Nb5d.CKO,
setname = "score.IPAN1_CKOup")## Summarizing data# scoring DEGs.IPAN2
GEX.seur <- add_geneset_score(obj = GEX.seur,
Assay = "SCT",
geneset = test1_up.list$IPAN2$Nb5d.PBS,
setname = "score.IPAN2_PBSup")## Summarizing dataGEX.seur <- add_geneset_score(obj = GEX.seur,
Assay = "SCT",
geneset = test1_up.list$IPAN2$Nb5d.INF,
setname = "score.IPAN2_INFup")## Summarizing dataGEX.seur <- add_geneset_score(obj = GEX.seur,
Assay = "SCT",
geneset = test2_up.list$IPAN2$Nb5d.CTL,
setname = "score.IPAN2_CTLup")## Summarizing dataGEX.seur <- add_geneset_score(obj = GEX.seur,
Assay = "SCT",
geneset = test2_up.list$IPAN2$Nb5d.CKO,
setname = "score.IPAN2_CKOup")## Summarizing data# All_PBSup
x = "All_PBSup"
VlnPlot(GEX.seur, features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.25,0.45)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.25, 0.25, 0.2),
size=2.75
)
# All_INFup
x = "All_INFup"
VlnPlot(GEX.seur, features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.25,0.45)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.25, 0.25, 0.2),
size=2.75
)
# All_CTLup
x = "All_CTLup"
VlnPlot(GEX.seur, features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.25,0.45)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.25, 0.25, 0.2),
size=2.75
)
# All_CKOup
x = "All_CKOup"
VlnPlot(GEX.seur, features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.25,0.45)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.25, 0.25, 0.2),
size=2.75
)
vln.DEGs_IPAN1 <- list()
# IPAN1_PBSup
x = "IPAN1_PBSup"
vln.DEGs_IPAN1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno1 %in% c("IPAN1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0.1,1.15)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(1, 1, 0.85),
size=2.75
)
# IPAN1_INFup
x = "IPAN1_INFup"
vln.DEGs_IPAN1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno1 %in% c("IPAN1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,0.75)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.55, 0.55, 0.6),
size=2.75
)
# IPAN1_CTLup
x = "IPAN1_CTLup"
vln.DEGs_IPAN1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno1 %in% c("IPAN1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,0.75)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.55, 0.55, 0.58),
size=2.75
)
# IPAN1_CKOup
x = "IPAN1_CKOup"
vln.DEGs_IPAN1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno1 %in% c("IPAN1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0.1,1.45)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(1.35, 1.35, 1.15),
size=2.75
)
vln.DEGs_IPAN1## $IPAN1_PBSup
##
## $IPAN1_INFup
##
## $IPAN1_CTLup
##
## $IPAN1_CKOup
vln.DEGs_IPAN1.1 <- list()
# IPAN1_PBSup
x = "IPAN1_PBSup"
vln.DEGs_IPAN1.1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN1.1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1.1 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0.1,1.15)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(1, 1, 0.85),
size=2.75
)
# IPAN1_INFup
x = "IPAN1_INFup"
vln.DEGs_IPAN1.1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN1.1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1.1 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,0.65)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.55, 0.55, 0.6),
size=2.75
)
# IPAN1_CTLup
x = "IPAN1_CTLup"
vln.DEGs_IPAN1.1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN1.1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1.1 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,0.65)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.55, 0.55, 0.58),
size=2.75
)
# IPAN1_CKOup
x = "IPAN1_CKOup"
vln.DEGs_IPAN1.1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN1.1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1.1 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0.1,1.5)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(1.35, 1.35, 1.15),
size=2.75
)
vln.DEGs_IPAN1.1## $IPAN1_PBSup
##
## $IPAN1_INFup
##
## $IPAN1_CTLup
##
## $IPAN1_CKOup
vln.DEGs_IPAN1.2 <- list()
# IPAN1_PBSup
x = "IPAN1_PBSup"
vln.DEGs_IPAN1.2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN1.2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1.2 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0.15,1.15)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(1, 1, 0.85),
size=2.75
)
# IPAN1_INFup
x = "IPAN1_INFup"
vln.DEGs_IPAN1.2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN1.2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1.2 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,0.65)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.55, 0.55, 0.6),
size=2.75
)
# IPAN1_CTLup
x = "IPAN1_CTLup"
vln.DEGs_IPAN1.2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN1.2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1.2 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,0.65)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.55, 0.55, 0.58),
size=2.75
)
# IPAN1_CKOup
x = "IPAN1_CKOup"
vln.DEGs_IPAN1.2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN1.2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN1.2 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0.2,1.5)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(1.35, 1.35, 1.15),
size=2.75
)
vln.DEGs_IPAN1.2## $IPAN1_PBSup
##
## $IPAN1_INFup
##
## $IPAN1_CTLup
##
## $IPAN1_CKOup
vln.DEGs_IPAN2 <- list()
# IPAN2_PBSup
x = "IPAN2_PBSup"
vln.DEGs_IPAN2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno1 %in% c("IPAN2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,1.25)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(1.1, 1, 0.85),
size=2.75
)
# IPAN2_INFup
x = "IPAN2_INFup"
vln.DEGs_IPAN2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno1 %in% c("IPAN2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,0.85)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.7, 0.65, 0.75),
size=2.75
)
# IPAN2_CTLup
x = "IPAN2_CTLup"
vln.DEGs_IPAN2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno1 %in% c("IPAN2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,0.85)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.72, 0.72, 0.78),
size=2.75
)
# IPAN2_CKOup
x = "IPAN2_CKOup"
vln.DEGs_IPAN2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno1 %in% c("IPAN2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & #labs(title = "") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,1.85)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(1.65, 1.65, 1.55),
size=2.75
)
vln.DEGs_IPAN2## $IPAN2_PBSup
##
## $IPAN2_INFup
##
## $IPAN2_CTLup
##
## $IPAN2_CKOup
vln.DEGs_IPAN2.1 <- list()
# IPAN2_PBSup
x = "IPAN2_PBSup"
vln.DEGs_IPAN2.1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN2.1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN2.1 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,1.25)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(1.1, 1, 0.85),
size=2.75
)
# IPAN2_INFup
x = "IPAN2_INFup"
vln.DEGs_IPAN2.1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN2.1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN2.1 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,0.85)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.7, 0.65, 0.75),
size=2.75
)
# IPAN2_CTLup
x = "IPAN2_CTLup"
vln.DEGs_IPAN2.1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN2.1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN2.1 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,0.85)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.72, 0.72, 0.78),
size=2.75
)
# IPAN2_CKOup
x = "IPAN2_CKOup"
vln.DEGs_IPAN2.1[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN2.1")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN2.1 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,1.85)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(1.65, 1.65, 1.55),
size=2.75
)
vln.DEGs_IPAN2.1## $IPAN2_PBSup
##
## $IPAN2_INFup
##
## $IPAN2_CTLup
##
## $IPAN2_CKOup
vln.DEGs_IPAN2.2 <- list()
# IPAN2_PBSup
x = "IPAN2_PBSup"
vln.DEGs_IPAN2.2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN2.2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN2.2 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,1.1)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.95, 0.85, 0.75),
size=2.75
)
# IPAN2_INFup
x = "IPAN2_INFup"
vln.DEGs_IPAN2.2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN2.2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN2.2 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,0.65)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.5, 0.54, 0.58),
size=2.75
)
# IPAN2_CTLup
x = "IPAN2_CTLup"
vln.DEGs_IPAN2.2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN2.2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN2.2 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,0.7)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.5, 0.45, 0.58),
size=2.75
)
# IPAN2_CKOup
x = "IPAN2_CKOup"
vln.DEGs_IPAN2.2[[x]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c("IPAN2.2")), features = paste0("score.",x), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "IPAN2.2 only") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.1,1.65)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(1.45, 1.45, 1.35),
size=2.75
)
vln.DEGs_IPAN2.2## $IPAN2_PBSup
##
## $IPAN2_INFup
##
## $IPAN2_CTLup
##
## $IPAN2_CKOup
check IEGs
IEGs.1 <- c("Fos","Jun","Egr1","Arc",
"Npas4","Homer1","Nptx2","Crem",
"Syne1", # CPG2
"Ptgs2", # COX-2
"Nrn1", # CPG15
#"Nptx2", # Narp
"Bdnf",
"Plk2", # SNK
"Plat" # tPA
)
IEGs.1 <- IEGs.1[IEGs.1 %in% rownames(GEX.seur)]
IEGs.1## [1] "Fos" "Jun" "Egr1" "Arc" "Npas4" "Homer1" "Crem" "Syne1"
## [9] "Nrn1" "Bdnf" "Plk2" "Plat"length(IEGs.1)## [1] 12sum(IEGs.1 %in% rownames(GEX.seur))## [1] 12pm.All_new.IEGs <- DotPlot(GEX.seur, features = IEGs.1, group.by = "intAnno2",
cols = c("midnightblue","darkorange1")) +
theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1,size = 9.15))+ scale_y_discrete(limits=rev) + labs(title="All - intAnno2")
pm.All_new.IEGs
# scoring IEGs
GEX.seur <- add_geneset_score(obj = GEX.seur,
Assay = "SCT",
geneset = IEGs.1,
setname = "score.IEGs")## Summarizing dataneur.clusters## [1] "EMN1" "EMN2" "EMN3" "EMN4" "EMN5" "IMN1" "IMN2"
## [8] "IMN3" "IMN4" "IN1" "IN2" "IN3" "IPAN1.1" "IPAN1.2"
## [15] "IPAN2.1" "IPAN2.2" "IPAN3" "IPAN4" VlnPlot(GEX.seur, features = paste0("score.","IEGs"), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = "All neurons") &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.15,0.3)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.15, 0.15, 0.2),
size=2.75
)## Warning: Removed 1 rows containing missing values (geom_point).
list_new## $All
## [1] "EMN1" "EMN2" "EMN3" "EMN4" "EMN5" "IMN1" "IMN2"
## [8] "IMN3" "IMN4" "IN1" "IN2" "IN3" "IPAN1.1" "IPAN1.2"
## [15] "IPAN2.1" "IPAN2.2" "IPAN3" "IPAN4"
##
## $EMNs
## [1] "EMN1" "EMN2" "EMN3" "EMN4" "EMN5"
##
## $IMNs
## [1] "IMN1" "IMN2" "IMN3" "IMN4"
##
## $IPAN1
## [1] "IPAN1.1" "IPAN1.2"
##
## $IPAN2
## [1] "IPAN2.1" "IPAN2.2"
##
## $EMN1
## [1] "EMN1"
##
## $EMN2
## [1] "EMN2"
##
## $EMN3
## [1] "EMN3"
##
## $EMN4
## [1] "EMN4"
##
## $EMN5
## [1] "EMN5"
##
## $IMN1
## [1] "IMN1"
##
## $IMN2
## [1] "IMN2"
##
## $IMN3
## [1] "IMN3"
##
## $IMN4
## [1] "IMN4"
##
## $IN1
## [1] "IN1"
##
## $IN2
## [1] "IN2"
##
## $IN3
## [1] "IN3"
##
## $IPAN1.1
## [1] "IPAN1.1"
##
## $IPAN1.2
## [1] "IPAN1.2"
##
## $IPAN2.1
## [1] "IPAN2.1"
##
## $IPAN2.2
## [1] "IPAN2.2"
##
## $IPAN3
## [1] "IPAN3"
##
## $IPAN4
## [1] "IPAN4"length(list_new)## [1] 23vln.IEGs <- list()
for(NN in names(list_new)){
vln.IEGs[[NN]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c(list_new[[NN]])), features = paste0("score.","IEGs"), group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = NN) &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(-0.15,0.3)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.15, 0.15, 0.2),
size=2.75
)
}
#vln.IEGscowplot::plot_grid(
plotlist = vln.IEGs,
ncol = 6
)
vln.Syne1 <- list()
for(NN in names(list_new)){
vln.Syne1[[NN]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c(list_new[[NN]])), features = "Syne1", group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = NN) &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,2.5)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(1.25, 1.25, 1.6),
size=2.75
)
}
#vln.Syne1cowplot::plot_grid(
plotlist = vln.Syne1,
ncol = 6
)
vln.Fos <- list()
for(NN in names(list_new)){
vln.Fos[[NN]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c(list_new[[NN]])), features = "Fos", group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = NN) &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,1.5)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.25, 0.25, 0.5),
size=2.75
)
}
#vln.Fosvln.Jun <- list()
for(NN in names(list_new)){
vln.Jun[[NN]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c(list_new[[NN]])), features = "Jun", group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = NN) &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,1.5)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.25, 0.25, 0.5),
size=2.75
)
}
#vln.Junvln.Nrn1 <- list()
for(NN in names(list_new)){
vln.Nrn1[[NN]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c(list_new[[NN]])), features = "Nrn1", group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = NN) &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,1.5)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.25, 0.25, 0.5),
size=2.75
)
}
#vln.Nrn1vln.Bdnf <- list()
for(NN in names(list_new)){
vln.Bdnf[[NN]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c(list_new[[NN]])), features = "Bdnf", group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = NN) &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,1.5)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.25, 0.25, 0.5),
size=2.75
)
}
#vln.Bdnfvln.Il13ra1 <- list()
for(NN in names(list_new)){
vln.Il13ra1[[NN]] <- VlnPlot(subset(GEX.seur, subset= intAnno2 %in% c(list_new[[NN]])), features = "Il13ra1", group.by = "cnt", cols = color.cnt, pt.size = 0) +
NoLegend() &
geom_jitter(alpha=0.1, shape=16, width = 0.2, size = 0.15) & labs(x = NN) &
geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & coord_cartesian(ylim = c(0,2)) &
ggpubr::stat_compare_means(aes(lable = ..p.signif..),
method = "wilcox.test",
comparisons = list(c("Nb5d.PBS","Nb5d.INF"),
c("Nb5d.CTL","Nb5d.CKO"),
c("Nb5d.INF","Nb5d.CTL")),
label.y = c(0.85, 0.85, 0.95),
size=2.75
)
}
#vln.Il13ra1#saveRDS(GEX.seur, "./sn10x_WYS.sct_anno.s.rds")
#GEX.seur <- readRDS("./sn10x_WYS.sct_anno.s.rds")